!= 主成分相変化
!
!= Phase change of atmospheric major component
!
! Authors::   Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: major_comp_phase_change.f90,v 1.3 2025/09/17 22:00:00 takepiro Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module major_comp_phase_change
  !
  != 主成分相変化
  !
  != Phase change of atmospheric major component
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  !== Procedures List
  ! 
!!$  ! DryConvAdjust :: 乾燥対流調節
!!$  ! ------------  :: ------------
!!$  ! DryConvAdjust :: Dry convective adjustment
  !
  !== NAMELIST
  !
  ! NAMELIST#major_comp_phase_change_nml
  !

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

  ! 格子点設定
  ! 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

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING     ! 文字列.       Strings. 

  ! NAMELIST ファイル入力に関するユーティリティ
  ! Utilities for NAMELIST file input
  !
  use namelist_util, only: MaxNmlArySize
                              ! NAMELIST から読み込む配列の最大サイズ. 
                              ! Maximum size of arrays loaded from NAMELIST

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

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private


  ! 公開手続き
  ! Public procedure
  !
  public :: MajorCompPhaseChangeInAtm
  public :: MajorCompPhaseChangeOnGround
  public :: MajorCompPhaseChangeInit


  ! 公開変数
  ! Public variables
  !
  logical, save, public:: major_comp_phase_change_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  ! 非公開変数
  ! Private variables
  !
  logical, save :: FlagMajCompPhaseChange
  logical, save :: FlagModMom


  character(*), parameter:: module_name = 'major_comp_phase_change'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: major_comp_phase_change.f90,v 1.3 2025/09/17 22:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version


contains

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

  subroutine MajorCompPhaseChangeInAtmTest(  &
    & xyr_Press, xyz_Press, xyz_Height,      &  ! (in)
    & xy_Ps, xyz_Temp, xy_SurfMajCompIce     &  ! (inout)
    & )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

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

    ! 物理定数設定
    ! 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

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only :     &
      & SaturateMajorCompCondTemp,      &
      & SaturateMajorCompInqLatentHeat


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(in   ):: xyz_Height(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! 
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyz_DelAtmMass      (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Atmospheric mass in a layer
    real(DP):: xy_FallingIce       (0:imax-1, 1:jmax)
                              !
                              !
    real(DP):: xyz_DelMajCompIce   (0:imax-1, 1:jmax, 1:kmax)
                              !
                              !
    real(DP):: xy_DelSurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              !
    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency

    real(DP):: xyz_TempCond   (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xy_SurfTempCond(0:imax-1, 1:jmax)
    real(DP):: SpecHeatCO2Ice

    real(DP):: LatentHeatMajCompSubl

    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

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

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


    if ( .not. FlagMajCompPhaseChange ) return


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


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! 調節前 "Temp" の保存
    ! Store "Temp" before adjustment
    !
    xyz_TempB    = xyz_Temp


    call SaturateMajorCompCondTemp(  &
      & xyz_Press,                   & ! (in)
      & xyz_TempCond                 & ! (inout)
      & )


    do k = 1, kmax
      xyz_DelAtmMass(:,:,k) = ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) / Grav
    end do


    k = kmax
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyz_Temp(i,j,k) < xyz_TempCond(i,j,k) ) then
          xyz_DelMajCompIce(i,j,k) =                                 &
            &   CpDry * ( xyz_TempCond(i,j,k) - xyz_Temp(i,j,k) )    &
            &   * xyz_DelAtmMass(i,j,k)                              &
            &   / LatentHeatMajCompSubl
          xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
        else
          xyz_DelMajCompIce(i,j,k) = 0.0_DP
        end if
      end do
    end do
    !
    xy_FallingIce = 0.0_DP
    do k = kmax-1, 1, -1
      xy_FallingIce = xy_FallingIce + xyz_DelMajCompIce(:,:,k+1)
      do j = 1, jmax
        do i = 0, imax-1
          SpecHeatCO2Ice = 349.0_DP + 4.8_DP * xyz_TempCond(i,j,k)
          !                                            Forget et al. (1998)
          xyz_DelMajCompIce(i,j,k) =                                     &
            &   CpDry * ( xyz_TempCond(i,j,k) - xyz_Temp(i,j,k) )        &
            &   * xyz_DelAtmMass(i,j,k)                                  &
            &   / LatentHeatMajCompSubl                                  &
            & - (   Grav * ( xyz_Height(i,j,k+1) - xyz_Height(i,j,k) )   &
            &     + SpecHeatCO2Ice                                       &
            &       * ( xyz_TempCond(i,j,k+1) - xyz_TempCond(i,j,k) ) )  &
            &   / LatentHeatMajCompSubl                                  &
            &   * xy_FallingIce(i,j)
          if ( ( xy_FallingIce(i,j) + xyz_DelMajCompIce(i,j,k) ) >= 0.0_DP ) then
            xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
          else
            xyz_DelMajCompIce(i,j,k) = - xy_FallingIce(i,j)
            xyz_Temp(i,j,k) = xyz_Temp(i,j,k)                                &
              & + ( - LatentHeatMajCompSubl                                  &
              &     + Grav * ( xyz_Height(i,j,k+1) - xyz_Height(i,j,k) )     &
              &     + SpecHeatCO2Ice                                         &
              &         * ( xyz_TempCond(i,j,k+1) - xyz_TempCond(i,j,k) ) )  &
              &     / ( CpDry * xyz_DelAtmMass(i,j,k) )                      &
              &     * xy_FallingIce(i,j)
          end if
        end do
      end do
    end do


    ! Ice falling on the surface
    !   This may result in supersaturation in the lowest level.
    !
    xy_FallingIce = xy_FallingIce + xyz_DelMajCompIce(:,:,1)
    k = 1
    do j = 1, jmax
      do i = 0, imax-1
        xy_DelSurfMajCompIce(i,j) =                                    &
          & - (   Grav * ( xyz_Height(i,j,1) - 0.0_DP ) )              &
          &   / LatentHeatMajCompSubl                                  &
          &   * xy_FallingIce(i,j)



          SpecHeatCO2Ice = 349.0_DP + 4.8_DP * xy_SurfTempCond(i,j)
          !                                            Forget et al. (1998)
          xy_DelSurfMajCompIce(i,j) =                                    &
!!$            &   CpDry * ( xyz_TempCond(i,j,k) - xyz_Temp(i,j,k) )        &
!!$            &   * xyz_DelAtmMass(i,j,k)                                  &
!!$            &   / LatentHeatMajCompSubl                                  &
            & - (   Grav * ( xyz_Height(i,j,1) - 0.0_DP )                &
            &     + SpecHeatCO2Ice                                       &
            &       * ( xyz_TempCond(i,j,1) - xy_SurfTempCond(i,j) ) )   &
            &   / LatentHeatMajCompSubl                                  &
            &   * xy_FallingIce(i,j)


        if ( ( xy_FallingIce(i,j) + xy_DelSurfMajCompIce(i,j) ) >= 0.0_DP ) then
          ! Part of ice sublimes.
          ! NOTE: In this case, temperature in the lowest layer should be 
          ! condensation temperature. So, actually, the set of temperature is 
          ! meaningless.
          xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
        else
          ! All falling ice sublimes.
          ! NOTE: The formulation below is different from that by Forget et al.
          ! (1998). The latent heat and heat by potential energy release and 
          ! heating ice is distributed in the lowest layer in this model, not
          ! to the soil.
          xy_DelSurfMajCompIce(i,j) = - xy_FallingIce(i,j)
          xyz_Temp(i,j,k) = xyz_Temp(i,j,k)                                &
            & + ( - LatentHeatMajCompSubl                                  &
            &     + Grav * ( xyz_Height(i,j,1) - 0.0_DP )                  &
            &     + SpecHeatCO2Ice                                         &
            &         * ( xyz_TempCond(i,j,1) - xy_SurfTempCond(i,j) ) )   &
            &     / ( CpDry * xyz_DelAtmMass(i,j,k) )                      &
            &     * xy_FallingIce(i,j)
       end if
       end do
    end do


    ! 温度変化率
    ! Calculate temperature tendency
    !
    xyz_DTempDt = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime )


    !
    ! Surface major component ice adjustment
    !
    xy_DSurfMajCompIceDt = 0.0_DP
    do k = kmax, 1, -1
      xy_DSurfMajCompIceDt = xy_DSurfMajCompIceDt + xyz_DelMajCompIce(:,:,k)
    end do
    xy_DSurfMajCompIceDt = xy_DSurfMajCompIceDt + xy_DelSurfMajCompIce(i,j)
    !
    xy_SurfMajCompIce = xy_SurfMajCompIce &
      & + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )


    !
    ! Surface pressure adjustment
    !
    xy_Ps = xy_Ps - xy_DSurfMajCompIceDt * Grav * ( 2.0_DP * DelTime )


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


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

  end subroutine MajorCompPhaseChangeInAtmTest

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

  subroutine MajorCompPhaseChangeInAtmBK(   &
    & xyr_Press, xyz_Press,               &  ! (in)
    & xy_Ps, xyz_Temp, xy_SurfMajCompIce  &  ! (inout)
    & )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

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

    ! 物理定数設定
    ! 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

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only :    &
      & SaturateMajorCompCondTemp,     &
      & SaturateMajorCompInqLatentHeat


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency

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

    real(DP):: LatentHeatMajCompSubl

    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

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

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


    if ( .not. FlagMajCompPhaseChange ) return


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


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! 調節前 "Temp" の保存
    ! Store "Temp" before adjustment
    !
    xyz_TempB    = xyz_Temp


    call SaturateMajorCompCondTemp( &
      & xyz_Press,                  & ! (in)
      & xyz_TempCond                & ! (inout)
      & )


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_Temp(i,j,k) < xyz_TempCond(i,j,k) ) then
            xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
          end if
        end do
      end do
    end do

    ! 温度変化率
    ! Calculate temperature tendency
    !
    xyz_DTempDt = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime )


    !
    ! Surface major component ice adjustment
    !
    xy_DSurfMajCompIceDt = 0.0_DP
    do k = kmax, 1, -1
      xy_DSurfMajCompIceDt = xy_DSurfMajCompIceDt                    &
        & + CpDry * xyz_DTempDt(:,:,k)                               &
        &   * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav       &
        &   / LatentHeatMajCompSubl
    end do
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )


    !
    ! Surface pressure adjustment
    !
    xy_Ps = xy_Ps - xy_DSurfMajCompIceDt * Grav * ( 2.0_DP * DelTime )


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


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

  end subroutine MajorCompPhaseChangeInAtmBK

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

  subroutine MajorCompPhaseChangeInAtmBK2(              &
    & xyr_Press, xyz_Press,                          & ! (in)
    & xy_Ps, xyz_Temp, xyzf_QMix, xy_SurfMajCompIce  & ! (inout)
    & )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

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

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : &
      & ncmax, &
      & IndexTKE

    ! 物理定数設定
    ! 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

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVars

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only :    &
      & SaturateMajorCompCondTemp,     &
      & SaturateMajorCompInqLatentHeat

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only: MassFixerColumn

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyzf_QMix        (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: LatentHeatMajCompSubl

    real(DP):: xy_PsB              (0:imax-1, 1:jmax)
    real(DP):: xy_PsA              (0:imax-1, 1:jmax)
    real(DP):: xyr_PressB          (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_PressA          (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyzf_QMixB          (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)


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

    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency
    real(DP):: xy_DPsDt            (0:imax-1, 1:jmax)

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

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

    integer :: mmax
    real(DP):: xyza_Array       (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    logical :: a_FlagSurfaceSink(1:ncmax)
    real(DP):: xyra_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax)

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

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

    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:: m
    integer:: n

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

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


    if ( .not. FlagMajCompPhaseChange ) return


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


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! Store variables
    !
    xyz_TempB  = xyz_Temp
    xyzf_QMixB = xyzf_QMix


    call SaturateMajorCompCondTemp( &
      & xyz_Press,                  & ! (in)
      & xyz_TempCond                & ! (inout)
      & )


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_Temp(i,j,k) < xyz_TempCond(i,j,k) ) then
            xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
          end if
        end do
      end do
    end do

    ! 温度変化率
    ! Calculate temperature tendency
    !
    xyz_DTempDt = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime )


    xyr_DPressDt(:,:,kmax) = 0.0_DP
    do k = kmax, 1, -1
      xy_DSurfMajCompIceDtOneLayer =                                 &
        &   CpDry * xyz_DTempDt(:,:,k)                               &
        &   * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav       &
        &   / LatentHeatMajCompSubl
      xyr_DPressDt(:,:,k-1) = xyr_DPressDt(:,:,k) &
        & - xy_DSurfMajCompIceDtOneLayer * Grav
    end do

    xy_DPsDt             = xyr_DPressDt(:,:,0)
    xy_DSurfMajCompIceDt = - xy_DPsDt / Grav



    ! packing
    mmax = ncmax
    do m = 1, mmax
      n = m
      xyza_Array(:,:,:,m) = xyzf_QMix(:,:,:,n)
      if ( n == IndexTKE ) then
        a_FlagSurfaceSink(m) = .true.
      else
        a_FlagSurfaceSink(m) = .false.
      end if
    end do

    call MajorCompPhaseChangeCalcFlow( &
      & xyr_Press, xyr_DPressDt,               & ! (in)
      & mmax, a_FlagSurfaceSink, xyza_Array,   & ! (in)
      & xyra_MassFlow                          & ! (out)
      & )

    ! unpacking
    do m = 1, mmax
      xyrf_MassFlow(:,:,:,m) = xyra_MassFlow(:,:,:,m)
    end do


    ! Adjustment
    !   preparation
    xy_PsB = xy_Ps
    xy_PsA = xy_PsB + xy_DPsDt * ( 2.0_DP * DelTime )

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    xyz_TempTmp    = 300.0_DP
    xyz_QH2OVapTmp =   0.0_DP
    call AuxVars( &
      & xy_PsB, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressB                  & ! (out) optional
      & )
    call AuxVars( &
      & xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressA                  & ! (out) optional
      & )
    do k = 1, kmax
      xyz_DelAtmMassB(:,:,k) = ( xyr_PressB(:,:,k-1) - xyr_PressB(:,:,k) ) / Grav
      xyz_DelAtmMassA(:,:,k) = ( xyr_PressA(:,:,k-1) - xyr_PressA(:,:,k) ) / Grav
    end do
    !   Atmospheric composition
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMix(:,:,k,n) =                                              &
          &   (   xyz_DelAtmMassB(:,:,k) * xyzf_QMix(:,:,k,n)             &
          &     - ( xyrf_MassFlow(:,:,k,n) - xyrf_MassFlow(:,:,k-1,n) ) ) &
          & / xyz_DelAtmMassA(:,:,k)
      end do
    end do


    ! Surface major component ice adjustment
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )
    ! Surface pressure adjustment
    xy_Ps = xy_PsA


    call AuxVars( &
      & xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressA                  & ! (out) optional
      & )
    ! 成分の質量の補正
    ! Fix masses of constituents
    !
    call MassFixerColumn( &
      & xyr_PressA, & ! (in)
      & xyzf_QMix   & ! (inout)
      & )

    ! Check
    call MajorCompPhaseChangeConsChk( &
      & a_FlagSurfaceSink,            & ! (in)
      & xyz_DelAtmMassB, xyzf_QMixB,  & ! (in)
      & xyz_DelAtmMassA, xyzf_QMix    & ! (in)
      & )

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


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

  end subroutine MajorCompPhaseChangeInAtmBK2

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

  subroutine MajorCompPhaseChangeInAtm(              &
    & xyr_Press, xyz_Press,                          & ! (in)
    & xy_Ps, xyz_Temp, xyzf_QMix, xyz_U, xyz_V,      & ! (inout)
    & xy_SurfMajCompIce                              & ! (inout)
    & )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

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

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : &
      & ncmax, &
      & IndexTKE

    ! 物理定数設定
    ! 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

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVars

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only :    &
      & SaturateMajorCompCondTemp,     &
      & SaturateMajorCompInqLatentHeat

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only: MassFixerColumn

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyzf_QMix        (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(inout):: xyz_U            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xyz_V            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: LatentHeatMajCompSubl

    real(DP):: xy_PsB              (0:imax-1, 1:jmax)
    real(DP):: xy_PsA              (0:imax-1, 1:jmax)
    real(DP):: xyr_PressB          (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_PressA          (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyzf_QMixB          (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)

    real(DP):: xyz_DelAtmMass      (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xy_TempTmp          (0:imax-1, 1:jmax)
    real(DP):: xy_DTempDtCond      (0:imax-1, 1:jmax)
    real(DP):: xyr_AtmMassFallFlux (0:imax-1, 1:jmax, 0:kmax)

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

    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency
    real(DP):: xy_DPsDt            (0:imax-1, 1:jmax)

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

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

    integer :: mmax
    real(DP):: xyza_Array       (0:imax-1, 1:jmax, 1:kmax, 1:ncmax+1+1)
    logical :: a_FlagSurfaceSink                          (1:ncmax+1+1)
    real(DP):: xyra_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax+1+1)

    real(DP):: xyrf_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
    real(DP):: xyr_MomXFlow     (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_MomYFlow     (0:imax-1, 1:jmax, 0:kmax)

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

    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:: m
    integer:: n

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

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


    if ( .not. FlagMajCompPhaseChange ) return


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


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! Store variables
    !
    xyz_TempB  = xyz_Temp
    xyzf_QMixB = xyzf_QMix


    call SaturateMajorCompCondTemp( &
      & xyz_Press,                  & ! (in)
      & xyz_TempCond                & ! (inout)
      & )


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

    xyr_AtmMassFallFlux(:,:,kmax) = 0.0_DP
    do k = kmax, 1, -1
      ! sublimation of falling condensate
!!$      xy_DTempDtSubl = &
!!$        & - LatentHeatMajCompSubl * xyr_AtmMassFallFlux(:,:,k) &
!!$        &   / ( CpDry * xyz_DelAtmMass(:,:,k) )
!!$      xyz_Temp(:,:,k) = xyz_Temp(:,:,k) + xy_DTempDtSubl * ( 2.0_DP * DelTime )
!!$      xyr_AtmMassFallFlux(:,:,k) = 0.0_DP

      ! condensation
      xy_TempTmp = xyz_Temp(:,:,k)
      xyz_Temp(:,:,k) = max( xyz_TempCond(:,:,k), xyz_Temp(:,:,k) )
      xy_DTempDtCond = ( xyz_Temp(:,:,k) - xy_TempTmp ) / ( 2.0_DP * DelTime )
      xy_DSurfMajCompIceDtOneLayer =                           &
        &   CpDry * xyz_DelAtmMass(:,:,k) * xy_DTempDtCond     &
        &   / LatentHeatMajCompSubl
      xyr_AtmMassFallFlux(:,:,k-1) = xyr_AtmMassFallFlux(:,:,k) &
        & + xy_DSurfMajCompIceDtOneLayer
    end do

    xyr_DPressDt = - xyr_AtmMassFallFlux * Grav


    xy_DPsDt             = xyr_DPressDt(:,:,0)
    xy_DSurfMajCompIceDt = - xy_DPsDt / Grav



    ! packing
    if ( FlagModMom ) then
      mmax = ncmax + 1 + 1
    else
      mmax = ncmax
    end if
    do m = 1, ncmax
      n = m
      xyza_Array(:,:,:,m) = xyzf_QMix(:,:,:,n)
      if ( n == IndexTKE ) then
        a_FlagSurfaceSink(m) = .true.
      else
        a_FlagSurfaceSink(m) = .false.
      end if
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_U
      a_FlagSurfaceSink(m) = .true.
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_V
      a_FlagSurfaceSink(m) = .true.
    end if

    call MajorCompPhaseChangeCalcFlow( &
      & xyr_Press, xyr_DPressDt,                                     & ! (in)
      & mmax, a_FlagSurfaceSink(1:mmax), xyza_Array(:,:,:,1:mmax),   & ! (in)
      & xyra_MassFlow(:,:,:,1:mmax)                                  & ! (out)
      & )

    ! unpacking
    do m = 1, ncmax
      xyrf_MassFlow(:,:,:,m) = xyra_MassFlow(:,:,:,m)
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyr_MomXFlow = xyra_MassFlow(:,:,:,m)
      m = m + 1
      xyr_MomYFlow = xyra_MassFlow(:,:,:,m)
    end if


    ! Adjustment
    !   preparation
    xy_PsB = xy_Ps
    xy_PsA = xy_PsB + xy_DPsDt * ( 2.0_DP * DelTime )

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    xyz_TempTmp    = 300.0_DP
    xyz_QH2OVapTmp =   0.0_DP
    call AuxVars( &
      & xy_PsB, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressB                  & ! (out) optional
      & )
    call AuxVars( &
      & xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressA                  & ! (out) optional
      & )
    do k = 1, kmax
      xyz_DelAtmMassB(:,:,k) = ( xyr_PressB(:,:,k-1) - xyr_PressB(:,:,k) ) / Grav
      xyz_DelAtmMassA(:,:,k) = ( xyr_PressA(:,:,k-1) - xyr_PressA(:,:,k) ) / Grav
    end do
    !   Atmospheric composition
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMix(:,:,k,n) =                                              &
          &   (   xyz_DelAtmMassB(:,:,k) * xyzf_QMix(:,:,k,n)             &
          &     - ( xyrf_MassFlow(:,:,k,n) - xyrf_MassFlow(:,:,k-1,n) ) ) &
          & / xyz_DelAtmMassA(:,:,k)
      end do
    end do
    if ( FlagModMom ) then
      do k = 1, kmax
        ! Zonal wind
        xyz_U(:,:,k) =                                              &
          &   (   xyz_DelAtmMassB(:,:,k) * xyz_U(:,:,k)             &
          &     - ( xyr_MomXFlow(:,:,k) - xyr_MomXFlow(:,:,k-1) ) ) &
          & / xyz_DelAtmMassA(:,:,k)
        ! Meridional wind
        xyz_V(:,:,k) =                                              &
          &   (   xyz_DelAtmMassB(:,:,k) * xyz_V(:,:,k)             &
          &     - ( xyr_MomYFlow(:,:,k) - xyr_MomYFlow(:,:,k-1) ) ) &
          & / xyz_DelAtmMassA(:,:,k)
      end do
    end if


    ! Surface major component ice adjustment
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )
    ! Surface pressure adjustment
    xy_Ps = xy_PsA


    call AuxVars( &
      & xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressA                  & ! (out) optional
      & )
    ! 成分の質量の補正
    ! Fix masses of constituents
    !
    call MassFixerColumn( &
      & xyr_PressA, & ! (in)
      & xyzf_QMix   & ! (inout)
      & )

    ! Check
    call MajorCompPhaseChangeConsChk( &
      & a_FlagSurfaceSink,            & ! (in)
      & xyz_DelAtmMassB, xyzf_QMixB,  & ! (in)
      & xyz_DelAtmMassA, xyzf_QMix    & ! (in)
      & )

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


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

  end subroutine MajorCompPhaseChangeInAtm

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

  subroutine MajorCompPhaseChangeOnGround( &
    & xy_DPsDt, xy_DSurfMajCompIceDt,      & ! (in)
    & xy_Ps, xyzf_QMix, xyz_U, xyz_V,      & ! (inout)
    & xy_SurfMajCompIce                    & ! (inout)
    & )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

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

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : &
      & ncmax, &
      & IndexTKE

    ! 物理定数設定
    ! 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

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVars

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only: MassFixerColumn


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xy_DPsDt            (0:imax-1, 1:jmax)
    real(DP), intent(in   ):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
    real(DP), intent(inout):: xyzf_QMix        (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(inout):: xyz_U            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xyz_V            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xy_SurfMajCompIce   (0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_PsB              (0:imax-1, 1:jmax)
    real(DP):: xy_PsA              (0:imax-1, 1:jmax)
    real(DP):: xyr_PressB          (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_PressA          (0:imax-1, 1:jmax, 0:kmax)

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

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

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

    integer :: mmax
    real(DP):: xyza_Array       (0:imax-1, 1:jmax, 1:kmax, 1:ncmax+1+1)
    logical :: a_FlagSurfaceSink                          (1:ncmax+1+1)
    real(DP):: xyra_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax+1+1)

    real(DP):: xyrf_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
    real(DP):: xyr_MomXFlow     (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_MomYFlow     (0:imax-1, 1:jmax, 0:kmax)

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

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


    ! 実行文 ; Executable statement
    !

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


    if ( .not. FlagMajCompPhaseChange ) return


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


    xy_PsB = xy_Ps
    xy_PsA = xy_PsB + xy_DPsDt * ( 2.0_DP * DelTime )

    xyzf_QMixB = xyzf_QMix

    xyz_TempTmp    = 300.0_DP
    xyz_QH2OVapTmp =   0.0_DP

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    call AuxVars( &
      & xy_PsB, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressB                  & ! (out) optional
      & )
    call AuxVars( &
      & xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressA                  & ! (out) optional
      & )


    xyr_DPressDt = ( xyr_PressA - xyr_PressB ) / ( 2.0_DP * Deltime )


    ! packing
    if ( FlagModMom ) then
      mmax = ncmax + 1 + 1
    else
      mmax = ncmax
    end if
    do m = 1, ncmax
      n = m
      xyza_Array(:,:,:,m) = xyzf_QMix(:,:,:,n)
      if ( n == IndexTKE ) then
        a_FlagSurfaceSink(m) = .true.
      else
        a_FlagSurfaceSink(m) = .false.
      end if
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_U
      a_FlagSurfaceSink(m) = .true.
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_V
      a_FlagSurfaceSink(m) = .true.
    end if

    call MajorCompPhaseChangeCalcFlow( &
      & xyr_PressB, xyr_DPressDt,                                    & ! (in)
      & mmax, a_FlagSurfaceSink(1:mmax), xyza_Array(:,:,:,1:mmax),   & ! (in)
      & xyra_MassFlow(:,:,:,1:mmax)                                  & ! (out)
      & )

    ! unpacking
    do m = 1, ncmax
      xyrf_MassFlow(:,:,:,m) = xyra_MassFlow(:,:,:,m)
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyr_MomXFlow = xyra_MassFlow(:,:,:,m)
      m = m + 1
      xyr_MomYFlow = xyra_MassFlow(:,:,:,m)
    end if


    ! Adjustment
    !   preparation
    do k = 1, kmax
      xyz_DelAtmMassB(:,:,k) = ( xyr_PressB(:,:,k-1) - xyr_PressB(:,:,k) ) / Grav
      xyz_DelAtmMassA(:,:,k) = ( xyr_PressA(:,:,k-1) - xyr_PressA(:,:,k) ) / Grav
    end do
    !   Atmospheric composition
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMix(:,:,k,n) =                                              &
          &   (   xyz_DelAtmMassB(:,:,k) * xyzf_QMix(:,:,k,n)             &
          &     - ( xyrf_MassFlow(:,:,k,n) - xyrf_MassFlow(:,:,k-1,n) ) ) &
          & / xyz_DelAtmMassA(:,:,k)
      end do
    end do
    if ( FlagModMom ) then
      do k = 1, kmax
        ! Zonal wind
        xyz_U(:,:,k) =                                              &
          &   (   xyz_DelAtmMassB(:,:,k) * xyz_U(:,:,k)             &
          &     - ( xyr_MomXFlow(:,:,k) - xyr_MomXFlow(:,:,k-1) ) ) &
          & / xyz_DelAtmMassA(:,:,k)
        ! Meridional wind
        xyz_V(:,:,k) =                                              &
          &   (   xyz_DelAtmMassB(:,:,k) * xyz_V(:,:,k)             &
          &     - ( xyr_MomYFlow(:,:,k) - xyr_MomYFlow(:,:,k-1) ) ) &
          & / xyz_DelAtmMassA(:,:,k)
      end do
    end if

    !   Surface major component ice
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )
    !   Surface pressure
    xy_Ps = xy_PsA


    call AuxVars( &
      & xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp,    & ! (in )
      & xyr_Press = xyr_PressA                  & ! (out) optional
      & )
    ! 成分の質量の補正
    ! Fix masses of constituents
    !
    call MassFixerColumn( &
      & xyr_PressA, & ! (in)
      & xyzf_QMix   & ! (inout)
      & )


    ! Check
    call MajorCompPhaseChangeConsChk( &
      & a_FlagSurfaceSink,            & ! (in)
      & xyz_DelAtmMassB, xyzf_QMixB,  & ! (in)
      & xyz_DelAtmMassA, xyzf_QMix    & ! (in)
      & )


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


  end subroutine MajorCompPhaseChangeOnGround

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

  subroutine MajorCompPhaseChangeCalcFlow( &
    & xyr_Press, xyr_DPressDt,             & ! (in)
    & mmax, a_FlagSurfaceSink, xyza_Array, & ! (in)
    & xyra_MassFlow                        & ! (out)
    & )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

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

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in ):: xyr_Press    (0:imax-1, 1:jmax, 0:kmax)
                              ! pressure
    real(DP), intent(in ):: xyr_DPressDt (0:imax-1, 1:jmax, 0:kmax)
    integer , intent(in ):: mmax
    logical , intent(in ):: a_FlagSurfaceSink(1:mmax)
    real(DP), intent(in ):: xyza_Array   (0:imax-1, 1:jmax, 1:kmax, 1:mmax)
    real(DP), intent(out):: xyra_MassFlow(0:imax-1, 1:jmax, 0:kmax, 1:mmax)

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyr_DPPress  (0:imax-1, 1:jmax, 0:kmax)
                              ! pressure at departure point
    real(DP):: DelAtmMass
    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:: k2              ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: m


    ! 実行文 ; Executable statement
    !

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


    if ( .not. FlagMajCompPhaseChange ) return


    xyr_DPPress = xyr_Press + xyr_DPressDt * ( 2.0_DP * DelTime )

    ! check
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyr_DPPress(i,j,k-1) < xyr_DPPress(i,j,k) ) then
            call MessageNotify( 'E', module_name, 'Order of departure points are inappropriate, P(k=%d)=%f < P(k=%d)=%f.', &
              & i = (/ k-1, k /), d = (/ xyr_DPPress(i,j,k-1), xyr_DPPress(i,j,k) /) )
          end if
        end do
      end do
    end do

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

!!$          if ( xyr_DPressDt(i,j,k) >= 0.0_DP ) then
          if ( xyr_DPPress(i,j,k) > xyr_Press(i,j,k) ) then

            sum_upward_mass_transport : do k2 = k, 1, -1
              if ( xyr_DPPress(i,j,k) > xyr_Press(i,j,k2-1) ) then
                DelAtmMass = ( xyr_Press(i,j,k2-1) - xyr_Press(i,j,k2) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) &
                    & + xyza_Array(i,j,k2,m) * DelAtmMass
                end do
              else
                DelAtmMass = ( xyr_DPPress(i,j,k) - xyr_Press(i,j,k2) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) &
                    & + xyza_Array(i,j,k2,m) * DelAtmMass
                end do
                exit sum_upward_mass_transport
              end if
            end do sum_upward_mass_transport

          else

            sum_downward_mass_transport : do k2 = k+1, kmax
              if ( xyr_DPPress(i,j,k) < xyr_Press(i,j,k2  ) ) then
                DelAtmMass = ( xyr_Press(i,j,k2-1) - xyr_Press(i,j,k2) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) &
                    & - xyza_Array(i,j,k2,m) * DelAtmMass
                end do
              else
                DelAtmMass = ( xyr_Press(i,j,k2-1) - xyr_DPPress(i,j,k) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) &
                    & - xyza_Array(i,j,k2,m) * DelAtmMass
                end do
                exit sum_downward_mass_transport
              end if
            end do sum_downward_mass_transport

          end if

        end do
      end do

    end do


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

          ! not surface sink
          if ( xyr_DPressDt(i,j,k) <= 0.0_DP ) then
            do m = 1, mmax
              if ( .not. a_FlagSurfaceSink(m) ) then
                xyra_MassFlow(i,j,k,m) = 0.0_DP
              end if
            end do
          end if

        end do
      end do

    end do



  end subroutine MajorCompPhaseChangeCalcFlow

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

  subroutine MajorCompPhaseChangeConsChk( &
    & a_FlagSurfaceSink,                  & ! (in)
    & xyz_DelAtmMassB, xyzf_QMixB,        & ! (in)
    & xyz_DelAtmMassA, xyzf_QMixA         & ! (in)
    & )

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : ncmax

    ! 物理定数設定
    ! 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

    logical , intent(in) :: a_FlagSurfaceSink(1:ncmax)
    real(DP), intent(in) :: xyz_DelAtmMassB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyzf_QMixB     (0:imax-1, 1:jmax, 1:kmax, ncmax)
    real(DP), intent(in) :: xyz_DelAtmMassA(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyzf_QMixA     (0:imax-1, 1:jmax, 1:kmax, ncmax)

    ! Local variables
    !
    real(DP) :: ValB
    real(DP) :: ValA
    real(DP) :: xyf_SumB(0:imax-1, 1:jmax, 1:ncmax)
    real(DP) :: xyf_SumA(0:imax-1, 1:jmax, 1:ncmax)
    real(DP) :: Ratio
    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: n


    xyf_SumB = 0.0_DP
    xyf_SumA = 0.0_DP
    do n = 1, ncmax
      do k = kmax, 1, -1
        xyf_SumB(:,:,n) = xyf_SumB(:,:,n) &
          & + xyz_DelAtmMassB(:,:,k) * xyzf_QMixB(:,:,k,n)
        xyf_SumA(:,:,n) = xyf_SumA(:,:,n) &
          & + xyz_DelAtmMassA(:,:,k) * xyzf_QMixA(:,:,k,n)
      end do
    end do
    do n = 1, ncmax
      if ( .not. a_FlagSurfaceSink(n) ) then
        do j = 1, jmax
          do i = 0, imax-1
            ValB = xyf_SumB(i,j,n)
            ValA = xyf_SumA(i,j,n)

            Ratio = ( ValA - ValB ) / ( ValA + 1.0d-100 )
            if ( abs( Ratio ) > 1.0d-10 ) then
              if ( ( ValB < 0.0_DP ) .and. ( abs( ValB ) < 1.0e-20_DP ) ) then
                ! Do nothing
              else
                call MessageNotify( 'M', module_name, 'Mass No. %d is not conserved, %f.', i = (/ n /), d = (/ Ratio /) )
              end if
            end if
          end do
        end do
      end if
    end do


  end subroutine MajorCompPhaseChangeConsChk

!!$  subroutine MajorCompPhaseChangeLimitTemp( &
!!$    & xyr_Press, xyz_Press,  &  ! (in)
!!$    & xy_SurfTemp, xyz_Temp  &  ! (inout)
!!$    & )
!!$    !
!!$    ! CO2 相変化
!!$    !
!!$    ! CO2 phase change
!!$    !
!!$
!!$    ! モジュール引用 ; USE statements
!!$    !
!!$
!!$    ! 時刻管理
!!$    ! Time control
!!$    !
!!$    use timeset, only: &
!!$      & DelTime, &            ! $ \Delta t $
!!$      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
!!$      & TimesetClockStart, TimesetClockStop
!!$
!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    use gtool_historyauto, only: HistoryAutoPut
!!$
!!$
!!$    ! 宣言文 ; Declaration statements
!!$    !
!!$    implicit none
!!$
!!$    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
!!$                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
!!$                              ! Air pressure (half level)
!!$    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! $ p $ . 気圧 (整数レベル). 
!!$                              ! Air pressure (full level)
!!$    real(DP), intent(inout):: xy_SurfTemp(0:imax-1, 1:jmax)
!!$                              ! $ T_s $ .   惑星表面温度. Surface temperature
!!$    real(DP), intent(inout):: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! $ T $ .     温度. Temperature
!!$
!!$    ! 作業変数
!!$    ! Work variables
!!$    !
!!$    real(DP):: xy_SurfTempB  (0:imax-1, 1:jmax)
!!$                              ! 調節前の惑星表面温度. 
!!$                              ! Surface temperature before adjustment
!!$    real(DP):: xyz_TempB     (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! 調節前の温度. 
!!$                              ! Temperature before adjustment
!!$    real(DP):: xy_DSurfTempDt(0:imax-1, 1:jmax)
!!$                              ! 惑星表面温度変化率. 
!!$                              ! Surface temperature tendency
!!$    real(DP):: xyz_DTempDt   (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! 温度変化率. 
!!$                              ! Temperature tendency
!!$
!!$    real(DP):: xy_SurfTempCond(0:imax-1, 1:jmax)
!!$    real(DP):: xyz_TempCond   (0:imax-1, 1:jmax, 1:kmax)
!!$
!!$    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
!!$
!!$    logical :: FlagCheckPs
!!$
!!$
!!$    ! 実行文 ; Executable statement
!!$    !
!!$
!!$    ! 計算時間計測開始
!!$    ! Start measurement of computation time
!!$    !
!!$    call TimesetClockStart( module_name )
!!$
!!$    ! 初期化
!!$    ! Initialization
!!$    !
!!$    if ( .not. major_comp_phase_change_inited ) call MajorCompPhaseChangeInit
!!$
!!$    if ( .not. FlagUse ) return
!!$
!!$
!!$    FlagCheckPs = .false.
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xyr_Press(i,j,0) > 1.0d4 ) then
!!$          FlagCheckPs = .true.
!!$        end if
!!$      end do
!!$    end do
!!$    if ( FlagCheckPs ) then
!!$      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
!!$    end if
!!$
!!$
!!$    ! 調節前 "Temp" の保存
!!$    ! Store "Temp" before adjustment
!!$    !
!!$    xy_SurfTempB = xy_SurfTemp
!!$    xyz_TempB    = xyz_Temp
!!$
!!$
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        xy_SurfTempCond(i,j) = &
!!$          & 149.2d0 + 6.48d0 * log( 0.135d0 * xyr_Press(i,j,0) * 1.0d-2 )
!!$      end do
!!$    end do
!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          xyz_TempCond(i,j,k) = &
!!$            & 149.2d0 + 6.48d0 * log( 0.135d0 * xyz_Press(i,j,k) * 1.0d-2 )
!!$        end do
!!$      end do
!!$    end do
!!$
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_SurfTemp(i,j) < xy_SurfTempCond(i,j) ) then
!!$          xy_SurfTemp(i,j) = xy_SurfTempCond(i,j)
!!$        end if
!!$      end do
!!$    end do
!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xyz_Temp(i,j,k) < xyz_TempCond(i,j,k) ) then
!!$            xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
!!$          end if
!!$        end do
!!$      end do
!!$    end do
!!$
!!$
!!$    ! 温度変化率
!!$    ! Calculate temperature tendency
!!$    !
!!$    xy_DSurfTempDt = ( xy_SurfTemp - xy_SurfTempB ) / ( 2.0_DP * DelTime )
!!$    xyz_DTempDt    = ( xyz_Temp    - xyz_TempB    ) / ( 2.0_DP * DelTime )
!!$
!!$
!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    call HistoryAutoPut( TimeN, 'DSurfTempDtCO2PhaseChange', xy_DSurfTempDt )
!!$    call HistoryAutoPut( TimeN, 'DTempDtCO2PhaseChange'    , xyz_DTempDt    )
!!$
!!$
!!$    ! 計算時間計測一時停止
!!$    ! Pause measurement of computation time
!!$    !
!!$    call TimesetClockStop( module_name )
!!$
!!$  end subroutine MajorCompPhaseChangeLimitTemp

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

  subroutine MajorCompPhaseChangeInit(             &
    & ArgFlagMajCompPhaseChange, CondMajCompName   & ! (in)
    & )
    !
    ! major_comp_phase_change モジュールの初期化を行います. 
    ! NAMELIST#major_comp_phase_change_nml の読み込みはこの手続きで行われます. 
    !
    ! "major_comp_phase_change" module is initialized. 
    ! "NAMELIST#major_comp_phase_change_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

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

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

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

    ! 補助的な変数を計算するサブルーチン・関数群
    ! Subroutines and functions for calculating auxiliary variables
    !
    use auxiliary, only : AuxVarsInit

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only : &
      & SaturateMajorCompInit

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only : MassFixerInit

    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage


    ! 宣言文 ; Declaration statements
    !
    implicit none

    logical     , intent(in) :: ArgFlagMajCompPhaseChange
    character(*), intent(in) :: CondMajCompName


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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /major_comp_phase_change_nml/ &
      & FlagModMom

          ! デフォルト値については初期化手続 "major_comp_phase_change#MajorCompPhaseChangeInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "major_comp_phase_change#MajorCompPhaseChangeInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( major_comp_phase_change_inited ) return


    FlagMajCompPhaseChange = ArgFlagMajCompPhaseChange


    ! デフォルト値の設定
    ! Default values settings
    !
    FlagModMom = .false.


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

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = major_comp_phase_change_nml )
    end if


    if ( FlagMajCompPhaseChange ) then
      ! 主成分相変化
      ! Phase change of atmospheric major component
      !
      call SaturateMajorCompInit(  &
        & CondMajCompName          & ! (in)
        & )
    end if

    ! 補助的な変数を計算するサブルーチン・関数群
    ! Subroutines and functions for calculating auxiliary variables
    !
    call AuxVarsInit

    ! 質量の補正
    ! Mass fixer
    !
    call MassFixerInit


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'DSurfTempDtMajCompPhaseChange', &
      & (/ 'lon ', 'lat ', 'time' /),                             &
      & 'heating by major component phase change', 'K s-1' )
    call HistoryAutoAddVariable( 'DTempDtMajCompPhaseChange',     &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /),                 &
      & 'heating by major component phase change', 'K s-1' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  FlagModMom = %b', l = (/ FlagModMom /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    major_comp_phase_change_inited = .true.

  end subroutine MajorCompPhaseChangeInit

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

end module major_comp_phase_change