!= 地下熱伝導モデルを用いた場合の陰解法による時間積分
!
!= Time integration by using implicit scheme in case using subsurface thermal diffusion model
!
! Authors::   Yasuhiro Morikawa, Yukiko Yamada, Yoshiyuki O. Takahashi
! Version::   $Id: phy_implicit_sdh_V3.f90,v 1.1 2015/01/29 12:05:01 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2009. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module phy_implicit_sdh_V3
  !
  != 地下熱伝導モデルを用いた場合の陰解法による時間積分
  !
  != Time integration by using implicit scheme in case using subsurface thermal diffusion model
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  !== Procedures List
  !
  ! PhyImplSDHTendency              :: 時間変化率の計算
  ! PhyImplSDHSetMethodFromMatthews :: SurfType から計算法インデクスの作成
  ! PhyImplSDHInit                  :: 初期化
  ! ------------------------------- :: ------------
  ! PhyImplSDHTendency              :: Calculate tendency
  ! PhyImplSDHSetMethodFromMatthews :: Set index for calculation method
  ! PhyImplSDHInit                  :: Initialization
  !
  !--
  !== NAMELIST
  !
  ! NAMELIST#phy_implicit_sdh_V3_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
    &                  kslmax  ! 地下の鉛直層数. 
                               ! Number of subsurface vertical level

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

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

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

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: PhyImplSDHV3Tendency
  public :: PhyImplSDHV3SetMethodMatthews
  public :: PhyImplSDHV3CorSOTempBySnowMelt
  public :: PhyImplSDHV3Init


  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  logical, save      :: FlagBucketModel
                              ! flag for use of bucket model
  logical, save      :: FlagSnow
                              ! flag for treating snow
  logical, save      :: FlagSlabOcean
                              ! flag for use of slab ocean
  logical, save      :: FlagMajCompPhaseChange
                              ! flag for use of slab ocean

  logical, save      :: FlagSublimation
                              !
                              ! flag for treating sublimation

  integer, parameter :: IndexSpcMajComp = 0
  integer, parameter :: IndexSpcH2O     = 1

  integer, save      :: NumMaxItr
  real(DP), save     :: TempItrCrit


  ! Index for calculation method
  !
  integer, parameter :: IndexPresTs         = 10 ! Prescribed surface temperature
  integer, parameter :: IndexSlabOcean      = 11 ! Slab ocean
  integer, parameter :: IndexSeaIce         = 12 ! Sea ice
  integer, parameter :: IndexLand           = 13 ! Land
  integer, parameter :: IndexLandWithPresTs = 14 ! Land with prescribed surface temperature

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

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


  real(DP), save:: SOHeatCapacity
                              ! Slab ocean heat capacity (J m-2 K-1)


contains

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

  subroutine PhyImplSDHV3SetMethodMatthews(    &
    & xy_SurfCond, xy_SurfType, xy_SeaIceConc, & ! (in)
    & xy_IndexCalcMethod                       & ! (out)
    & )
    !
    !
    !
    ! Set index for calculation method from Matthews' index
    !

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

    ! 雪, 氷の割合
    ! snow/ice fraction
    !
    use snowice_frac, only : SeaIceAboveThreshold


    ! 宣言文 ; Declaration statements
    !

    integer , intent(in ) :: xy_SurfCond       (0:imax-1, 1:jmax)
                              ! 
                              ! Surface condition
    integer , intent(in ) :: xy_SurfType       (0:imax-1, 1:jmax)
                              ! 土地利用.
                              ! Surface index
    real(DP), intent(in ) :: xy_SeaIceConc     (0:imax-1, 1:jmax)
                              ! 海氷密度 (0 <= xy_SeaIceConc <= 1)
                              ! Sea ice concentration (0 <= xy_SeaIceConc <= 1)
    integer , intent(out) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method


    ! 作業変数
    ! Work variables
    !

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


    ! 実行文 ; Executable statement
    !

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


!!$    if ( .not. FlagSSModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$    end if

    ! FlagBucketModel は関係ないよね?
    ! SSModel 強制にした時点で, 水蒸気は地面と分離したから. 
!!$    if ( .not. FlagBucketModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$    end if


    !
    ! Set index for calculation method
    !
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_SurfCond(i,j) == 0 ) then
          if ( xy_SurfType(i,j) >= 1 ) then
            ! land
            xy_IndexCalcMethod(i,j) = IndexLandWithPresTs
          else
            ! prescribed surface temperature
            xy_IndexCalcMethod(i,j) = IndexPresTs
          end if
        else
          if ( xy_SurfType(i,j) >= 1 ) then
            ! land
            xy_IndexCalcMethod(i,j) = IndexLand
          else
            if ( SeaIceAboveThreshold( xy_SeaIceConc(i,j) ) ) then
              ! sea ice
              xy_IndexCalcMethod(i,j) = IndexSeaIce
            else if ( FlagSlabOcean ) then
              ! slab ocean
              xy_IndexCalcMethod(i,j) = IndexSlabOcean
            else
              ! prescribed surface temperature
!!$              xy_IndexCalcMethod(i,j) = IndexPresTs
              call MessageNotify( 'E', module_name, 'Unexpected behavior.' )
            end if
          end if
        end if
      end do
    end do


  end subroutine PhyImplSDHV3SetMethodMatthews

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

  subroutine PhyImplSDHV3CorSOTempBySnowMelt(  &
    & xy_IndexCalcMethod,                    & ! (in   )
    & xy_SurfSnowFlux,                       & ! (in   )
    & xy_SurfTemp                            & ! (inout)
    & )
    !
    !
    !
    ! Set index for calculation method from Matthews' index
    !

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & LatentHeatFusion
                              ! $ L $ [J kg-1] . 
                              ! 融解の潜熱. 
                              ! Latent heat of fusion

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


    ! 宣言文 ; Declaration statements
    !
    integer , intent(in   ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method
    real(DP), intent(in   ) :: xy_SurfSnowFlux   (0:imax-1, 1:jmax)
    real(DP), intent(inout) :: xy_SurfTemp       (0:imax-1, 1:jmax)


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_SeaIceThickness(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


    ! 実行文 ; Executable statement
    !

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


!!$    if ( .not. FlagSSModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$    end if

    ! FlagBucketModel は関係ないよね?
    ! SSModel 強制にした時点で, 水蒸気は地面と分離したから. 
!!$    if ( .not. FlagBucketModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$    end if


    !
    ! Set sea ice thickness
    !
    xy_SeaIceThickness = SeaIceThickness

    !
    ! Set index for calculation method
    !
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
        case ( IndexSeaIce )
          xy_SurfTemp(i,j) = xy_SurfTemp(i,j)                              &
            & + LatentHeatFusion * xy_SurfSnowFlux(i,j) * 2.0_DP * DelTime &
            &   / ( SeaIceVolHeatCap * xy_SeaIceThickness(i,j) )
        case ( IndexSlabOcean )
          xy_SurfTemp(i,j) = xy_SurfTemp(i,j)                              &
            & + LatentHeatFusion * xy_SurfSnowFlux(i,j) * 2.0_DP * DelTime &
            &   / SOHeatCapacity
        case ( IndexPresTs )
        case default
          call MessageNotify( 'E', module_name, 'This index is inappropriate.' )
        end select
      end do
    end do


  end subroutine PhyImplSDHV3CorSOTempBySnowMelt

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

  subroutine PhyImplSDHV3Tendency(                             &
    & xy_IndexCalcMethod, xy_BucketFlagOceanGrid,              & ! (in)
    & xy_SnowFrac,                                             & ! (in)
    & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in)
    & xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux,                & ! (out)
    & xyr_SoilHeatFlux,                                        & ! (in)
    & xyr_RadSFlux, xyr_RadLFlux,                              & ! (in)
    & xy_DeepSubSurfHeatFlux,                                  & ! (in)
    & xyz_TempB, xy_SurfTemp, xyz_SoilTemp,                    & ! (in)
    & xyzf_QMixB,                                              & ! (in)
    & xy_SurfHumidCoef,                                        & ! (in)
    & xy_SoilHeatCap, xy_SoilHeatDiffCoef,                     & ! (in)
    & xyra_DelRadLFlux,                                        & ! (in)
    & xyr_Press, xyz_Exner, xyr_Exner,                         & ! (in)
    & xyr_VirTemp, xyz_Height,                                 & ! (in)
    & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,     & ! (in)
    & xy_SurfVelTransCoef, xy_SurfTempTransCoef,               & ! (in)
    & xy_SurfQVapTransCoef,                                    & ! (in)
    & xyr_SoilTempTransCoef,                                   & ! (in)
    & xy_SurfMajCompIceB,                                      & ! (in)
    & xy_SoilMoistB, xy_SurfSnowB,                             & ! (in)
    & xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt,           & ! (out)
    & xy_DSurfTempDt,                                          & ! (out)
    & xyz_DSoilTempDt,                                         & ! (out)
    & xy_DPsDt, xy_DSurfMajCompIceDt,                          & ! (out)
    & xy_DSoilMoistDt,                                         & ! (out)
    & xy_DSurfSnowDt                                           & ! (out)
    & )
    !
    ! 時間変化率の計算を行います. 
    !
    ! Calculate tendencies. 
    !

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

    ! MPI 関連ルーチン
    ! MPI related routines
    !
    use mpi_wrapper, only: myrank, MPIWrapperFindMaxVal

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_SSDepth,   &    ! subsurface grid on interface of layer
      & z_SSDepth         ! subsurface grid at midpoint of layer

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

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: &
      & TempCondWater,       &
      & SnowVolHeatCap,      &
      & SnowDens,            &
      & SnowMaxThermDepth,   &
      & SeaIceVolHeatCap   , &
      & SeaIceThermCondCoef, &
      & SeaIceThickness,    &
      & TempBelowSeaIce

    ! 飽和比湿の算出
    ! Evaluation of saturation specific humidity
    !
    use saturate, only: &
      & xy_CalcQVapSatOnLiq, &
      & xy_CalcQVapSatOnSol, &
      & xy_CalcDQVapSatDTempOnLiq, &
      & xy_CalcDQVapSatDTempOnSol

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

    ! バケツモデル
    ! Bucket model
    !
    use Bucket_Model, only :            &
      & BucketModEvapAndLatentHeatFlux

    ! 地表面フラックスユーティリティ
    ! Surface flux utility routines
    !
    use surface_flux_util, only : SurfaceFluxUtilLimitFlux

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: xy_IndexCalcMethod(0:imax-1, 1:jmax) 
                              ! 
                              ! Index for calculation method
    logical , intent(in):: xy_BucketFlagOceanGrid(0:imax-1, 1:jmax)
                              !
                              ! Flag for ocean grid point used in bucket model
    real(DP), intent(in):: xy_SnowFrac  (0:imax-1, 1:jmax)
                              !
                              ! Snow fraction
    real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux

    real(DP), intent(out):: xy_SurfH2OVapFlux(0:imax-1, 1:jmax)
                              ! 惑星表面水蒸気フラックス.
                              ! Water vapor flux at the surface
    real(DP), intent(out):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface

    real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
                              ! 土壌中の熱フラックス (W m-2)
                              ! Heat flux in sub-surface soil (W m-2)

    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux

    real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
                              ! 地中熱フラックス. 
                              ! "Deep subsurface heat flux"
                              ! Heat flux at the bottom of surface/soil layer.
    real(DP), intent(in):: xyz_TempB(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度. 
                              ! Temperature
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xyz_SoilTemp (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度 (K)
                              ! Soil temperature (K)
    real(DP), intent(in):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! 
                              ! 
    real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP), intent(in ):: xy_SoilHeatCap (0:imax-1, 1:jmax)
                              ! 土壌熱容量 (J K-1 kg-1)
                              ! Specific heat of soil (J K-1 kg-1)
    real(DP), intent(in ):: xy_SoilHeatDiffCoef (0:imax-1, 1:jmax)
                              ! 土壌熱伝導係数 (J m-3 K-1)
                              ! Heat conduction coefficient of soil (J m-3 K-1)

    real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)
    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)

    real(DP), intent(in):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity

    real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
                              ! 輸送係数：土壌温度.
                              ! Transfer coefficient: soil temperature

    real(DP), intent(in):: xy_SurfMajCompIceB  (0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount.

    real(DP), intent(in):: xy_SoilMoistB (0:imax-1, 1:jmax)
                              ! 土壌水分.
                              ! Soil moisture.
    real(DP), intent(in):: xy_SurfSnowB (0:imax-1, 1:jmax)
                              ! 積雪量.
                              ! Surface snow amount.

    real(DP), intent(out):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{u}{t} $ . 東西風速変化. 
                              ! Eastward wind tendency
    real(DP), intent(out):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{v}{t} $ . 南北風速変化. 
                              ! Northward wind tendency
    real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(out):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 質量混合比変化. 
                              ! Mass mixing ratio tendency
    real(DP), intent(out):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)
    real(DP), intent(out):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
                              ! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
                              ! Temperature tendency (K s-1)

    real(DP), intent(out):: xy_DPsDt            (0:imax-1, 1:jmax)
    real(DP), intent(out):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)

    real(DP), intent(out):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
                              ! 土壌温度時間変化率 (kg m-2 s-1)
                              ! Soil temperature tendency (kg m-2 s-1)
    real(DP), intent(out):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
                              ! 積雪率時間変化率 (kg m-2 s-1)
                              ! Surface snow amount tendency (kg m-2 s-1)

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
                              ! 地表熱容量. 
                              ! Surface heat capacity

    real(DP) :: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              !
                              ! Sea ice thickness

    real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP) :: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass of constituents

    real(DP):: xyza_UVMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 速度陰解行列. 
                              ! Implicit matrix about velocity 
    real(DP):: xyz_UVec (0:imax-1, 1:jmax, 1:kmax)
                              ! 速度陰解ベクトル. 
                              ! Implicit vector about velocity 
    real(DP):: xyz_VVec (0:imax-1, 1:jmax, 1:kmax)
                              ! 速度陰解ベクトル. 
                              ! Implicit vector about velocity 
    real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP):: xyza_QMixMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 質量混合比陰解行列. 
                              ! Implicit matrix about mass mixing ratio
    real(DP):: xyzf_QMixVec(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! 質量混合比陰解ベクトル. 
                              ! Implicit vector about mass mixing ratio
    real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance

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

    real(DP):: xyza_UVLUMtx (0:imax-1, 1:jmax, 1:kmax,-1:1)
                              ! LU 行列. 
                              ! LU matrix

!!$    real(DP):: xyza_TempQVapLUMtx (0:imax-1, 1:jmax, -kmax:kmax, -1:1)
!!$                              ! LU 行列. 
!!$                              ! LU matrix
!!$    real(DP):: xyz_DelTempQVap (0:imax-1, 1:jmax, -kmax:kmax)
!!$                              ! $ T q $ の時間変化. 
!!$                              ! Tendency of $ T q $ 
!!$
!!$    real(DP):: xyza_TempLUMtx (0:imax-1, 1:jmax, 0:kmax, -1:1)
!!$                              ! LU 行列.
!!$                              ! LU matrix
!!$    real(DP):: xyz_DelTempLUVec (0:imax-1, 1:jmax, 0:kmax)
!!$                              ! $ T q $ の時間変化.
!!$                              ! Tendency of $ T q $
    real(DP):: xyza_QMixLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xyz_DelQMixLUVec (0:imax-1, 1:jmax, 1:kmax)
                              ! $ q $ の時間変化.
                              ! Tendency of $ q $

    real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature
    real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
                              ! $ T, Tg $ の時間変化.
                              ! Tendency of $ T $ and $ Tg |

    real(DP):: SurfSnowATentative
                              ! 積雪量の仮の値 (kg m-2)
                              ! pseudo value of surface snow amount (kg m-2)

    real(DP):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by major component ice sublimation
                              ! (variable only for debug)
    real(DP):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by melt
                              ! (variable only for debug)
    real(DP):: xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by sea ice melt
                              ! (variable only for debug)

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

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

    real(DP):: xy_SurfQVapSatOnLiq(0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface

    real(DP):: xy_SurfSoilHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面土壌熱伝導フラックス.
                              ! Soil heat conduction flux at the surface


    real(DP):: xyz_TempSave   (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_TempA      (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: MaxTempInc
    real(DP):: xy_SurfTempSave(0:imax-1, 1:jmax)
    real(DP):: xy_SurfTempA   (0:imax-1, 1:jmax)
    real(DP):: MaxSurfTempInc

    real(DP):: a_LocalMax (2)
    real(DP):: a_GlobalMax(2)

    integer:: iitr

    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               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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


    ! Check kskmax
    !
    if ( kslmax < 1 ) then
      call MessageNotify( 'E', module_name, 'kslmax is less than 1.' )
    end if

!!$    if ( .not. FlagSSModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$    end if

    ! FlagBucketModel は関係ないよね?
    ! SSModel 強制にした時点で, 水蒸気は地面と分離したから. 
!!$    if ( .not. FlagBucketModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$    end if

    ! Set heat capacity, tentatively here.
    !
!!$    xy_SurfHeatCapacity = 0.0_DP

    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          !   J K-1 m-3 kg m-2 (kg m-3)-1 = J K-1 m-5 m3 = J K-1 m-2
          xy_SurfHeatCapacity(i,j) = &
            & SnowVolHeatCap &
            & * min( max( xy_SurfSnowB(i,j) / SnowDens, 0.0_DP ), SnowMaxThermDepth )
        case default
          xy_SurfHeatCapacity(i,j) = 0.0_DP
        end select
      end do
    end do

    !
    ! Set sea ice thickness
    !
    xy_SeaIceThickness = SeaIceThickness


    ! Set coefficients for soil heat flux
    !
!!$    xyr_SoilTempTransCoef = xyr_BareSoilTempTransCoef
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_SurfSnowB(i,j) > 
!!$    xyr_SoilTempTransCoef = xyr_BareSoilTempTransCoef


    ! 輸送係数の計算
    ! Calculate transfer coefficient
    !
    xyr_VelTransCoef (:,:,0)    = 0.0_DP
    xyr_VelTransCoef (:,:,kmax) = 0.0_DP
    xyr_TempTransCoef(:,:,0)    = 0.0_DP
    xyr_TempTransCoef(:,:,kmax) = 0.0_DP
    xyr_QMixTransCoef(:,:,0)    = 0.0_DP
    xyr_QMixTransCoef(:,:,kmax) = 0.0_DP

    do k = 1, kmax-1
      xyr_VelTransCoef(:,:,k) =                                     &
        &   xyr_VelDiffCoef(:,:,k)                                  &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )

      xyr_TempTransCoef(:,:,k) =                                    &
        &   xyr_TempDiffCoef(:,:,k)                                 &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )

      xyr_QMixTransCoef(:,:,k) =                                    &
        &   xyr_QMixDiffCoef(:,:,k)                                 &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
    end do


    ! Calculation for momentum diffusion
    !
    call PhyImplSDHV3TendencyMomCore( &
      & xyr_MomFluxX, xyr_MomFluxY,         & ! (in)
      & xyr_Press,                          & ! (in)
      & xyr_VelTransCoef,                   & ! (in)
      & xy_SurfVelTransCoef,                & ! (in)
      & xyz_DUDt, xyz_DVDt                  & ! (out)
      & )


    ! Calculation for thermal diffusion
    !

    ! 飽和比湿の計算
    ! Calculate saturated specific humidity
    !
    xy_SurfQVapSatOnLiq       = &
      & xy_CalcQVapSatOnLiq      ( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSatOnSol       = &
      & xy_CalcQVapSatOnSol      ( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSat       = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq &
      & + xy_SnowFrac              * xy_SurfQVapSatOnSol
    xy_SurfDQVapSatDTempOnLiq = &
      & xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
    xy_SurfDQVapSatDTempOnSol = &
      & xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
    xy_SurfDQVapSatDTemp = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq &
      & + xy_SnowFrac              * xy_SurfDQVapSatDTempOnSol


    ! Initialization
    !
    xyzf_DQMixDt    = 0.0_DP
    xy_DSurfTempDt  = 0.0_DP
    xyz_DSoilTempDt = 0.0_DP

    xyz_TempSave    = xyz_TempB
    xy_SurfTempSave = xy_SurfTemp

    ! iteration
    iitr = 1
    loop_itr : do

      ! Tendencies of atmospheric and surface temperatures, and atmospheric 
      ! water vapor are solved with a fixed surface heat conduction flux at the surface. 
      ! Obtained tendencies of surface temperature and atmospheric water vapor will be 
      ! used below to estimate surface water vapor flux. 
      !
      xy_SurfSoilHeatFlux = xyr_SoilHeatFlux(:,:,0)                              &
        & - xyr_SoilTempTransCoef(:,:,0)                                         &
        &   * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0_DP * DelTime )
      !
      call PhyImplSDHV3TendencyHeatTQCore(                   &
        & xy_IndexCalcMethod,                                      & ! (in)
        & xy_SeaIceThickness,                                      & ! (in)
        & xy_SnowFrac,                                             & ! (in)
        & xyr_HeatFlux, xyrf_QMixFlux,                             & ! (in)
        & xy_SurfSoilHeatFlux,                                     & ! (in)
        & xyr_RadSFlux, xyr_RadLFlux,                              & ! (in)
        & xy_DeepSubSurfHeatFlux,                                  & ! (in)
        & xy_SurfTemp,                                             & ! (in)
        & xy_SurfHumidCoef,                                        & ! (in)
        & xy_SurfHeatCapacity,                                     & ! (in)
        & xyra_DelRadLFlux,                                        & ! (in)
        & xyr_Press, xyz_Exner, xyr_Exner,                         & ! (in)
        & xyr_VelTransCoef, xyr_TempTransCoef,                     & ! (in)
        & xyr_QMixTransCoef,                                       & ! (in)
        & xy_SurfVelTransCoef, xy_SurfTempTransCoef,               & ! (in)
        & xy_SurfQVapTransCoef,                                    & ! (in)
        & xyz_DTempDt, xyzf_DQMixDt,                               & ! (out)
        & xy_DSurfTempDt                                           & ! (out)
        & )


      ! Tendencies of atmospheric, surface, and subsurface temperatures, and 
      ! atmospheric water vapor are solved with a fixed surface water vapor flux. 
      !
      n = IndexH2OVap
      xy_SurfH2OVapFlux = xyrf_QMixFlux(:,:,0,n)               &
        & - xy_SurfHumidCoef * xy_SurfQVapTransCoef            &
        &   * ( xyzf_DQMixDt(:,:,1,n)                          &
        &     - xy_SurfDQVapSatDTemp * xy_DSurfTempDt )        &
        &   * 2.0_DP * DelTime
      ! Limit surface flux not to be negative atmospheric content
      ! IMPORTANT : Now, only the water vapor flux is restricted.
      call SurfaceFluxUtilLimitFlux( &
        & ( 2.0_DP * DelTime ),            & ! (in)
        & xyzf_QMixB, xyr_Press,           & ! (in)
        & xy_SurfH2OVapFlux &!, xyrf_QMixFlux(:,:,0,n)  & ! (inout)
        & )
      ! Calculation of latent heat flux
      xy_SurfLatentHeatFlux = LatentHeat * xy_SurfH2OVapFlux
      !
      if ( FlagBucketModel ) then
        ! バケツモデルのための地表面フラックス修正
        ! Modification of surface flux for bucket model
        !
        call BucketModEvapAndLatentHeatFlux(                         &
          & xy_BucketFlagOceanGrid, xy_SoilMoistB, xy_SurfSnowB,     & ! (in   )
          & xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux                 & ! (inout)
          & )
      end if
      !
      call PhyImplSDHV3TendencyHeatCore(                     &
        & xy_IndexCalcMethod,                                      & ! (in)
        & xy_SeaIceThickness,                                      & ! (in)
        & xyr_HeatFlux, xyrf_QMixFlux,                             & ! (in)
        & xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux,                & ! (in)
        & xyr_SoilHeatFlux,                                        & ! (in)
        & xyr_RadSFlux, xyr_RadLFlux,                              & ! (in)
        & xy_DeepSubSurfHeatFlux,                                  & ! (in)
        & xy_SurfTemp, xyz_SoilTemp,                               & ! (in)
        & xy_SurfHumidCoef,                                        & ! (in)
        & xy_SurfHeatCapacity,                                     & ! (in)
        & xy_SoilHeatCap, xy_SoilHeatDiffCoef,                     & ! (in)
        & xyra_DelRadLFlux,                                        & ! (in)
        & xyr_Press, xyz_Exner, xyr_Exner,                         & ! (in)
        & xyr_VelTransCoef, xyr_TempTransCoef,                     & ! (in)
        & xyr_QMixTransCoef,                                       & ! (in)
        & xy_SurfVelTransCoef, xy_SurfTempTransCoef,               & ! (in)
        & xy_SurfQVapTransCoef,                                    & ! (in)
        & xyr_SoilTempTransCoef,                                   & ! (in)
        & xy_SurfMajCompIceB,                                      & ! (in)
        & xy_SoilMoistB, xy_SurfSnowB,                             & ! (in)
        & xyz_DTempDt, xyzf_DQMixDt,                               & ! (out)
        & xy_DSurfTempDt,                                          & ! (out)
        & xyz_DSoilTempDt,                                         & ! (out)
        & xy_DPsDt, xy_DSurfMajCompIceDt,                          & ! (out)
        & xy_DSoilMoistDt,                                         & ! (out)
        & xy_DSurfSnowDt                                           & ! (out)
        & )


      ! Check

      xyz_TempA    = xyz_TempB   + xyz_DTempDt    * 2.0_DP * DelTime
      xy_SurfTempA = xy_SurfTemp + xy_DSurfTempDt *          DelTime

      MaxTempInc     = maxval( abs( xyz_TempA    - xyz_TempSave    ) )
      MaxSurfTempInc = maxval( abs( xy_SurfTempA - xy_SurfTempSave ) )

      a_LocalMax(1) = MaxTempInc
      a_LocalMax(2) = MaxSurfTempInc
      call MPIWrapperFindMaxVal( &
        & 2, a_LocalMax,         & ! (in)
        & a_GlobalMax            & ! (out)
        & )
      MaxTempInc     = a_GlobalMax(1)
      MaxSurfTempInc = a_GlobalMax(2)

!!$      if ( myrank == 0 ) then
!!$        call MessageNotify( 'M', module_name, &
!!$          & 'Itr: %d : dT = %f, dTs = %f',    &
!!$          & i = (/iitr/), d = (/ MaxTempInc, MaxSurfTempInc /) )
!!$      end if

      if ( ( MaxTempInc <= TempItrCrit ) .and. ( MaxSurfTempInc <= TempItrCrit ) ) then
        exit loop_itr
      end if

      xyz_TempSave    = xyz_TempA
      xy_SurfTempSave = xy_SurfTempA

      iitr = iitr + 1
      if ( iitr > NumMaxItr ) then
        if ( NumMaxItr > 2 ) then
          if ( myrank == 0 ) then
            call MessageNotify( 'M', module_name, &
              & 'Too many iterations, Itr: %d : dT = %f, dTs = %f', &
              & i = (/iitr/), d = (/ MaxTempInc, MaxSurfTempInc /) )
          end if
        end if
        exit loop_itr
      end if
    end do loop_itr

!!$    if ( myrank == 0 ) then
!!$      call MessageNotify( 'M', module_name, &
!!$        & 'Itr: %d : dT = %f, dTs = %f', &
!!$        & i = (/iitr/), d = (/ MaxTempInc, MaxSurfTempInc /) )
!!$    end if


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


  end subroutine PhyImplSDHV3Tendency

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

  subroutine PhyImplSDHV3TendencyMomCore( &
    & xyr_MomFluxX, xyr_MomFluxY,         & ! (in)
    & xyr_Press,                          & ! (in)
    & xyr_VelTransCoef,                   & ! (in)
    & xy_SurfVelTransCoef,                & ! (in)
    & xyz_DUDt, xyz_DVDt                  & ! (out)
    & )
    !
    ! 時間変化率の計算を行います. 
    !
    ! Calculate tendencies. 
    !

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

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

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

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux

    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)

    real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity

    real(DP), intent(out):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{u}{t} $ . 東西風速変化. 
                              ! Eastward wind tendency
    real(DP), intent(out):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{v}{t} $ . 南北風速変化. 
                              ! Northward wind tendency


    ! 作業変数
    ! Work variables
    !
    real(DP):: xyza_UVMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 速度陰解行列. 
                              ! Implicit matrix about velocity 
    real(DP):: xyz_UVec (0:imax-1, 1:jmax, 1:kmax)
                              ! 速度陰解ベクトル. 
                              ! Implicit vector about velocity 
    real(DP):: xyz_VVec (0:imax-1, 1:jmax, 1:kmax)
                              ! 速度陰解ベクトル. 
                              ! Implicit vector about velocity 

    real(DP):: xyza_UVLUMtx (0:imax-1, 1:jmax, 1:kmax,-1:1)
                              ! LU 行列. 
                              ! LU matrix


    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               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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



!!$    if ( .not. FlagSSModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$    end if

    ! FlagBucketModel は関係ないよね?
    ! SSModel 強制にした時点で, 水蒸気は地面と分離したから. 
!!$    if ( .not. FlagBucketModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$    end if



    ! 陰解法のための行列作成
    ! Create matrices for implicit scheme
    !

    ! 鉛直拡散スキームの輸送係数から陰解行列の計算 (速度)
    ! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (velocity)
    !
    k = 1
    xyza_UVMtx  (:,:,k,-1) = 0.0_DP
    xyza_UVMtx  (:,:,k, 0) =                                                  &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xy_SurfVelTransCoef(:,:)                                            &
      & + xyr_VelTransCoef(:,:,k  )
    xyza_UVMtx  (:,:,k, 1) = &
      & - xyr_VelTransCoef(:,:,k)

    do k = 2, kmax-1
      xyza_UVMtx  (:,:,k,-1) = &
        & - xyr_VelTransCoef(:,:,k-1)
      xyza_UVMtx  (:,:,k, 0) =                                                  &
        & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + xyr_VelTransCoef(:,:,k-1)                                           &
        & + xyr_VelTransCoef(:,:,k  )
      xyza_UVMtx  (:,:,k, 1) = &
        & - xyr_VelTransCoef(:,:,k)
    end do

    k = kmax
    xyza_UVMtx  (:,:,k,-1) = &
      & - xyr_VelTransCoef(:,:,k-1)
    xyza_UVMtx  (:,:,k, 0) =                                                  &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xyr_VelTransCoef(:,:,k-1)
    xyza_UVMtx  (:,:,k, 1) = 0.0_DP

    do k = 1, kmax
      xyz_UVec(:,:,k) = - ( xyr_MomFluxX(:,:,k) - xyr_MomFluxX(:,:,k-1) )
      xyz_VVec(:,:,k) = - ( xyr_MomFluxY(:,:,k) - xyr_MomFluxY(:,:,k-1) )
    end do

    ! 東西風速, 南北風速の計算
    ! Calculate eastward and northward wind
    !
    xyza_UVLUMtx = xyza_UVMtx

    call PhyImplLUDecomp3( &
      & xyza_UVLUMtx, &     ! (inout)
      & imax * jmax, kmax ) ! (in)

    do k = 1, kmax
      xyz_DUDt(:,:,k) = xyz_UVec(:,:,k)
      xyz_DVDt(:,:,k) = xyz_VVec(:,:,k)
    end do

    call PhyImplLUSolve3( &
      & xyz_DUDt, &            ! (inout)
      & xyza_UVLUMtx, &        ! (in)
      & 1, imax * jmax, kmax ) ! (in)

    call PhyImplLUSolve3( &
      & xyz_DVDt, &            ! (inout)
      & xyza_UVLUMtx, &        ! (in)
      & 1, imax * jmax, kmax ) ! (in)

    do k = 1, kmax
      xyz_DUDt(:,:,k) = xyz_DUDt(:,:,k) / ( 2. * DelTime )
      xyz_DVDt(:,:,k) = xyz_DVDt(:,:,k) / ( 2. * DelTime )
    end do


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

  end subroutine PhyImplSDHV3TendencyMomCore

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

  subroutine PhyImplSDHV3TendencyHeatCore(                     &
    & xy_IndexCalcMethod,                                      & ! (in)
    & xy_SeaIceThickness,                                      & ! (in)
    & xyr_HeatFlux, xyrf_QMixFlux,                             & ! (in)
    & xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux,                & ! (in)
    & xyr_SoilHeatFlux,                                        & ! (in)
    & xyr_RadSFlux, xyr_RadLFlux,                              & ! (in)
    & xy_DeepSubSurfHeatFlux,                                  & ! (in)
    & xy_SurfTemp, xyz_SoilTemp,                               & ! (in)
    & xy_SurfHumidCoef,                                        & ! (in)
    & xy_SurfHeatCapacity,                                     & ! (in)
    & xy_SoilHeatCap, xy_SoilHeatDiffCoef,                     & ! (in)
    & xyra_DelRadLFlux,                                        & ! (in)
    & xyr_Press, xyz_Exner, xyr_Exner,                         & ! (in)
    & xyr_VelTransCoef, xyr_TempTransCoef,                     & ! (in)
    & xyr_QMixTransCoef,                                       & ! (in)
    & xy_SurfVelTransCoef, xy_SurfTempTransCoef,               & ! (in)
    & xy_SurfQVapTransCoef,                                    & ! (in)
    & xyr_SoilTempTransCoef,                                   & ! (in)
    & xy_SurfMajCompIceB,                                      & ! (in)
    & xy_SoilMoistB, xy_SurfSnowB,                             & ! (in)
    & xyz_DTempDt, xyzf_DQMixDt,                               & ! (out)
    & xy_DSurfTempDt,                                          & ! (out)
    & xyz_DSoilTempDt,                                         & ! (out)
    & xy_DPsDt, xy_DSurfMajCompIceDt,                          & ! (out)
    & xy_DSoilMoistDt,                                         & ! (out)
    & xy_DSurfSnowDt                                           & ! (out)
    & )
    !
    ! 時間変化率の計算を行います. 
    !
    ! Calculate tendencies. 
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_SSDepth,   &    ! subsurface grid on interface of layer
      & z_SSDepth         ! subsurface grid at midpoint of layer

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

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

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

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method

    real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              !
                              ! Sea ice thickness

    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux

    real(DP), intent(in):: xy_SurfH2OVapFlux(0:imax-1, 1:jmax)
                              ! 惑星表面水蒸気フラックス.
                              ! Water vapor flux at the surface
    real(DP), intent(in):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface

    real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
                              ! 土壌中の熱フラックス (W m-2)
                              ! Heat flux in sub-surface soil (W m-2)

    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux

    real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
                              ! 地中熱フラックス. 
                              ! "Deep subsurface heat flux"
                              ! Heat flux at the bottom of surface/soil layer.
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xyz_SoilTemp (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度 (K)
                              ! Soil temperature (K)
    real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
                              ! 地表熱容量. 
                              ! Surface heat capacity
    real(DP), intent(in ):: xy_SoilHeatCap (0:imax-1, 1:jmax)
                              ! 土壌熱容量 (J K-1 kg-1)
                              ! Specific heat of soil (J K-1 kg-1)
    real(DP), intent(in ):: xy_SoilHeatDiffCoef (0:imax-1, 1:jmax)
                              ! 土壌熱伝導係数 (J m-3 K-1)
                              ! Heat conduction coefficient of soil (J m-3 K-1)

    real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(in):: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass of constituents

    real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity

    real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
                              ! 輸送係数：土壌温度.
                              ! Transfer coefficient: soil temperature

    real(DP), intent(in):: xy_SurfMajCompIceB  (0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount.

    real(DP), intent(in):: xy_SoilMoistB(0:imax-1, 1:jmax)
    real(DP), intent(in):: xy_SurfSnowB (0:imax-1, 1:jmax)
                              ! 積雪量.
                              ! Surface snow amount.

    real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(out):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 質量混合比変化. 
                              ! Mass mixing ratio tendency
    real(DP), intent(out):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)
    real(DP), intent(out):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
                              ! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
                              ! Temperature tendency (K s-1)

    real(DP), intent(out):: xy_DPsDt            (0:imax-1, 1:jmax)
    real(DP), intent(out):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)

    real(DP), intent(out):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
                              ! 土壌温度時間変化率 (kg m-2 s-1)
                              ! Soil temperature tendency (kg m-2 s-1)
    real(DP), intent(out):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
                              ! 積雪率時間変化率 (kg m-2 s-1)
                              ! Surface snow amount tendency (kg m-2 s-1)

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP):: xyza_QMixMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 質量混合比陰解行列. 
                              ! Implicit matrix about mass mixing ratio
    real(DP):: xyzf_QMixVec(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! 質量混合比陰解ベクトル. 
                              ! Implicit vector about mass mixing ratio
    real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance

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


!!$    real(DP):: xyza_TempQVapLUMtx (0:imax-1, 1:jmax, -kmax:kmax, -1:1)
!!$                              ! LU 行列. 
!!$                              ! LU matrix
!!$    real(DP):: xyz_DelTempQVap (0:imax-1, 1:jmax, -kmax:kmax)
!!$                              ! $ T q $ の時間変化. 
!!$                              ! Tendency of $ T q $ 
!!$
!!$    real(DP):: xyza_TempLUMtx (0:imax-1, 1:jmax, 0:kmax, -1:1)
!!$                              ! LU 行列.
!!$                              ! LU matrix
!!$    real(DP):: xyz_DelTempLUVec (0:imax-1, 1:jmax, 0:kmax)
!!$                              ! $ T q $ の時間変化.
!!$                              ! Tendency of $ T q $
    real(DP):: xyza_QMixLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xyz_DelQMixLUVec (0:imax-1, 1:jmax, 1:kmax)
                              ! $ q $ の時間変化.
                              ! Tendency of $ q $

    real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature
    real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
                              ! $ T, Tg $ の時間変化.
                              ! Tendency of $ T $ and $ Tg |

    real(DP):: SurfSnowATentative
                              ! 積雪量の仮の値 (kg m-2)
                              ! pseudo value of surface snow amount (kg m-2)

    real(DP):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by major component ice sublimation
                              ! (variable only for debug)
    real(DP):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by melt
                              ! (variable only for debug)
    real(DP):: xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by sea ice melt
                              ! (variable only for debug)

    real(DP):: xy_SurfMajCompLiqB      (0:imax-1, 1:jmax)
    real(DP):: xy_LatHeatFluxByOtherSpc(0:imax-1, 1:jmax)

    real(DP):: xy_DAtmMassDt(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               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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



!!$    if ( .not. FlagSSModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$    end if

    ! FlagBucketModel は関係ないよね?
    ! SSModel 強制にした時点で, 水蒸気は地面と分離したから. 
!!$    if ( .not. FlagBucketModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$    end if



    ! 陰解法のための行列作成
    ! Create matrices for implicit scheme
    !

    ! 鉛直拡散スキームの輸送係数から陰解行列の計算 (温度)
    ! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (temperature)
    !
    k = 1
    xyza_TempMtx(:,:,k,-1) = &
      & - CpDry * xy_SurfTempTransCoef(:,:)
    xyza_TempMtx(:,:,k, 0) =                                                          &
      & - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k  ) * xy_SurfTempTransCoef(:,:) &
      & + CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k  )
    xyza_TempMtx(:,:,k, 1) = &
      & - CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k  )

    do k = 2, kmax-1
      xyza_TempMtx(:,:,k,-1) = &
        & - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
      xyza_TempMtx(:,:,k, 0) =                                                          &
        & - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k-1)&
        & + CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k  )
      xyza_TempMtx(:,:,k, 1) = &
        & - CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k  )
    end do

    k = kmax
    xyza_TempMtx(:,:,k,-1) = &
      & - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
    xyza_TempMtx(:,:,k, 0) =                                                          &
      & - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k-1)
    xyza_TempMtx(:,:,k, 1) = 0.0_DP

    do k = 1, kmax
      xyz_TempVec(:,:,k) = - ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) )
    end do


    ! 鉛直拡散スキームの輸送係数から陰解行列の計算 (比湿)
    ! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (specific humidity)
    !

    k = 1
    xyza_QMixMtx(:,:,k,-1) =                                                  &
      & 0.0_DP
    xyza_QMixMtx(:,:,k, 0) =                                                  &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xyr_QMixTransCoef(:,:,k  )
    xyza_QMixMtx(:,:,k, 1) =                                                  &
      & - xyr_QMixTransCoef(:,:,k  )

    do k = 2, kmax-1
      xyza_QMixMtx(:,:,k,-1) =                                                  &
        & - xyr_QMixTransCoef(:,:,k-1)
      xyza_QMixMtx(:,:,k, 0) =                                                  &
        & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + xyr_QMixTransCoef(:,:,k-1)                                          &
        & + xyr_QMixTransCoef(:,:,k  )
      xyza_QMixMtx(:,:,k, 1) =                                                  &
        & - xyr_QMixTransCoef(:,:,k  )
    end do

    k = kmax
    xyza_QMixMtx(:,:,k,-1) =                                                  &
      & - xyr_QMixTransCoef(:,:,k-1)
    xyza_QMixMtx(:,:,k, 0) =                                                  &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xyr_QMixTransCoef(:,:,k-1)
    xyza_QMixMtx(:,:,k, 1) = 0.0_DP

    do n = 1, ncmax
      if ( n == IndexH2OVap ) then
        do k = 1, 1
          xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xy_SurfH2OVapFlux )
        end do
        do k = 1+1, kmax
          xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
        end do
      else
        do k = 1, kmax
          xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
        end do
      end if
    end do



    ! 土壌温度計算用の輸送係数から陰解行列の計算 (土壌温度)
    ! Calculate implicit matrices by using transfer coefficient (soil temperature)
    !
    if ( kslmax /= 0 ) then ! xyr_SoilTempMtx is not used when kslmax = 0.

      do k = 1, kslmax-1
        xyaa_SoilTempMtx(:,:,k,-1) =                                                  &
          & - xyr_SoilTempTransCoef(:,:,k-1)
        xyaa_SoilTempMtx(:,:,k, 0) =                                                  &
          &   xy_SoilHeatCap(:,:) * ( r_SSDepth(k) - r_SSDepth(k-1) )                 &
          &     / ( 2.0_DP * DelTime )&
          & + xyr_SoilTempTransCoef(:,:,k-1)                                          &
          & + xyr_SoilTempTransCoef(:,:,k  )
        xyaa_SoilTempMtx(:,:,k, 1) =                                                  &
          & - xyr_SoilTempTransCoef(:,:,k  )
      end do

      k = kslmax
      xyaa_SoilTempMtx(:,:,k,-1) =                                   &
        & - xyr_SoilTempTransCoef(:,:,k-1)
      xyaa_SoilTempMtx(:,:,k, 0) =                                   &
        &   xy_SoilheatCap(:,:) * ( r_SSDepth(k) - r_SSDepth(k-1) )  &
        &     / ( 2.0_DP * DelTime )                                     &
        & + xyr_SoilTempTransCoef(:,:,k-1)
      xyaa_SoilTempMtx(:,:,k, 1) = 0.0_DP

    end if

    do k = 1, kslmax
      xya_SoilTempVec (:,:,k) = - ( xyr_SoilHeatFlux(:,:,k) - xyr_SoilHeatFlux(:,:,k-1) )
    end do


    ! 地表面過程の輸送係数から陰解行列の計算
    ! Calculate implicit matrices from transfer coefficient of surface process
    !
    do i = 0, imax-1
      do j = 1, jmax
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          xyaa_SurfMtx(i,j,0,-1) =                                             &
            &   xyr_SoilTempTransCoef(i,j,0)
          xyaa_SurfMtx(i,j,0, 0) =                                             &
            &   xy_SurfHeatCapacity(i,j) / ( 2.0_DP * DelTime )                    &
            & + CpDry * xy_SurfTempTransCoef(i,j)                              &
            & + xyra_DelRadLFlux(i,j,0,0)                                      &
            & - xyr_SoilTempTransCoef(i,j,0)
          xyaa_SurfMtx(i,j,0, 1) =                                             &
            & - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1)                    &
            &   * xy_SurfTempTransCoef(i,j)                                    &
            & + xyra_DelRadLFlux(i,j,0,1)
        case ( IndexSeaIce )
          ! sea ice
          xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
          xyaa_SurfMtx(i,j,0, 0) =                                        &
            &   SeaIceVolHeatCap * xy_SeaIceThickness(i,j)                &
            &     / ( 2.0_DP * DelTime ) &
            & + CpDry * xy_SurfTempTransCoef(i,j)                         &
            & + xyra_DelRadLFlux(i,j,0,0)                                 &
            & + SeaIceThermCondCoef / xy_SeaIceThickness(i,j)
          xyaa_SurfMtx(i,j,0, 1) = &
            & - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1)                    &
            &   * xy_SurfTempTransCoef(i,j)                                    &
            & + xyra_DelRadLFlux(i,j,0,1)
        case ( IndexSlabOcean )
          ! slab ocean ocean
          xyaa_SurfMtx(i,j,0,-1) =                                             &
            & 0.0_DP
          xyaa_SurfMtx(i,j,0, 0) =                                             &
            &   SOHeatCapacity / ( 2.0_DP * DelTime )                          &
            & + CpDry * xy_SurfTempTransCoef(i,j)                              &
            & + xyra_DelRadLFlux(i,j,0,0)
          xyaa_SurfMtx(i,j,0, 1) =                                             &
            & - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1)                    &
            &   * xy_SurfTempTransCoef(i,j)                                    &
            & + xyra_DelRadLFlux(i,j,0,1)
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
          xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
          xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select

      end do
    end do

    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          xy_SurfRH(i,j) =                                    &
            & - xyr_RadSFlux(i,j,0)                           &
            & - xyr_RadLFlux(i,j,0)                           &
            & - xyr_HeatFlux(i,j,0)                           &
            & - xy_SurfLatentHeatFlux(i,j)                    &
!            & + xy_DeepSubSurfHeatFlux(i,j)
            & + xyr_SoilHeatFlux(i,j,0)
        case ( IndexSeaIce )
          ! sea ice
          xy_SurfRH(i,j) =                                      &
            & - xyr_RadSFlux(i,j,0)                             &
            & - xyr_RadLFlux(i,j,0)                             &
            & - xyr_HeatFlux(i,j,0)                             &
            & - xy_SurfLatentHeatFlux(i,j)                      &
!              & + xy_DeepSubSurfHeatFlux(i,j)
!              & + xyr_SoilHeatFlux(i,j,0)
            & - SeaIceThermCondCoef                             &
            &   * ( xy_SurfTemp(i,j) - TempBelowSeaIce ) / xy_SeaIceThickness(i,j)
        case ( IndexSlabOcean )
          ! slab ocean
          xy_SurfRH(i,j) =                                       &
            & - xyr_RadSFlux(i,j,0)                              &
            & - xyr_RadLFlux(i,j,0)                              &
            & - xyr_HeatFlux(i,j,0)                              &
            & - xy_SurfLatentHeatFlux(i,j)                       !&
!              & + xy_DeepSubSurfHeatFlux(i,j)
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xy_SurfRH(i,j) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do


    ! 温度と比湿の計算
    ! Calculate temperature and specific humidity
    !

    do l = -1, 1
      do k = 1, kslmax
        xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
      end do
      k = 0
      xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
      do k = 1, kmax
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
      end do
    end do

    call PhyImplLUDecomp3(              &
      & xyaa_TempSoilTempLUMtx,         & ! (inout)
      & imax * jmax, kmax + 1 + kslmax  & ! (in)
      )

    do k = 1, kslmax
      xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
    end do
    k = 0
    xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
    do k = 1, kmax
      xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
    end do

    call PhyImplLUSolve3(                   &
      & xya_DelTempSoilTempLUVec,           & ! (inout)
      & xyaa_TempSoilTempLUMtx,             & ! (in)
      & 1, imax * jmax , kmax + 1 + kslmax  & ! (in)
      )

    do k = 1, kslmax
      do j = 1, jmax
        do i = 0, imax-1
          select case ( xy_IndexCalcMethod(i,j) )
          case ( IndexLand )
            xyz_DSoilTempDt(i,j,k) = &
              & xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
          case default
            xyz_DSoilTempDt(i,j,k) = 0.0_DP
          end select
        end do
      end do
    end do
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
        case ( IndexSeaIce )
          ! sea ice
          xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
        case ( IndexSlabOcean )
          ! slab ocean
          xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xy_DSurfTempDt(i,j) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    do k = 1, kmax
      xyz_DTempDt(:,:,k) = xya_DelTempSoilTempLUVec(:,:,k) / ( 2.0_DP * DelTime )
    end do


    !
    ! Calculation of tendencies of soil moisture and surface snow amount
    !
    if ( FlagBucketModel ) then
      if ( FlagSnow ) then
        ! Evaporation is subtracted from surface snow and soil moisture
        !
        do j = 1, jmax
          do i = 0, imax-1

!!$            if ( xyrf_QMixFlux(i,j,0,IndexH2OVap) >= 0.0_DP ) then
            if ( xy_SurfH2OVapFlux(i,j) >= 0.0_DP ) then

!!$              xy_DSurfSnowDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
              xy_DSurfSnowDt(i,j) = - xy_SurfH2OVapFlux(i,j)
              SurfSnowATentative = xy_SurfSnowB(i,j) &
                & + xy_DSurfSnowDt(i,j) * 2.0_DP * DelTime
              if ( SurfSnowATentative < 0.0_DP ) then
                xy_DSoilMoistDt(i,j) = SurfSnowATentative / ( 2.0_DP * DelTime )
                xy_DSurfSnowDt (i,j) = - xy_SurfSnowB(i,j) / ( 2.0_DP * DelTime )
              else
                xy_DSoilMoistDt(i,j) = 0.0_DP
              end if

            else

              if ( xy_SurfSnowB(i,j) > 0.0_DP ) then
!!$                xy_DSurfSnowDt (i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
                xy_DSurfSnowDt (i,j) = - xy_SurfH2OVapFlux(i,j)
                xy_DSoilMoistDt(i,j) = 0.0_DP
              else
                xy_DSurfSnowDt (i,j) = 0.0_DP
!!$                xy_DSoilMoistDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
                xy_DSoilMoistDt(i,j) = - xy_SurfH2OVapFlux(i,j)
              end if

            end if

          end do
        end do
      else
        ! Evaporation is subtracted from soil moisture
        !
!!$        xy_DSoilMoistDt = - xyrf_QMixFlux(:,:,0,IndexH2OVap)
        xy_DSoilMoistDt = - xy_SurfH2OVapFlux
        xy_DSurfSnowDt  = 0.0_DP
      end if
    else
      xy_DSoilMoistDt = 0.0_DP
      xy_DSurfSnowDt  = 0.0_DP
    end if


    ! Temporarily set
    !
    xy_DSurfMajCompIceDt = 0.0_DP


    if ( FlagMajCompPhaseChange ) then

      xy_DAtmMassDt        = 0.0_DP
      xy_DSurfMajCompIceDt = 0.0_DP
      ! Dummy values
      !
      xy_SurfMajCompLiqB       = 0.0_DP
      xy_LatHeatFluxByOtherSpc = 0.0_DP

      call PhyImplSDHV3IceSnowPhaseChgCor(               &
        & IndexSpcMajComp,                               & ! (in)
        & xyr_Press(:,:,0),                              & ! (in)
        & xyr_HeatFlux, xy_SurfLatentHeatFlux,           & ! (in)
        & xyr_SoilHeatFlux,                              & ! (in)
        & xyr_SoilTempTransCoef,                         & ! (in)
        & xyr_RadSFlux, xyr_RadLFlux,                    & ! (in)
        & xy_DeepSubSurfHeatFlux,                        & ! (in)
        & xy_SurfTemp, xyz_SoilTemp,                     & ! (in)
!        & xy_SurfSnowB,                                  & ! (in)
!        & xy_SurfLiqB, xy_SurfSolB,                      & ! (in)
        & xy_SurfMajCompLiqB, xy_SurfMajCompIceB,        & ! (in)
        & xy_SurfHeatCapacity,                           & ! (in)
        & xy_SoilHeatCap, xy_SoilHeatDiffCoef,           & ! (in)
        & xy_IndexCalcMethod,                            & ! (in)
        & xyra_DelRadLFlux,                              & ! (in)
        & xyz_Exner, xyr_Exner,                          & ! (in)
        & xy_SurfTempTransCoef,                          & ! (in)
        & xy_LatHeatFluxByOtherSpc,                      & ! (in)
        & xyza_TempMtx, xyz_TempVec,                     & ! (in)
        & xyaa_SurfMtx, xy_SurfRH,                       & ! (in)
        & xyaa_SoilTempMtx, xya_SoilTempVec,             & ! (in)
        & xyz_DTempDt,                                   & ! (in)
        & xy_DSurfTempDt,                                & ! (inout)
        & xyz_DSoilTempDt,                               & ! (inout)
!        & xy_DSoilMoistDt, xy_DSurfSnowDt,               & ! (inout)
!        & xy_LatHeatFluxBySnowMelt                       & ! (out)
        & xy_DAtmMassDt, xy_DSurfMajCompIceDt,           & ! (inout)
        & xy_LatHeatFluxByMajCompIceSubl                 & ! (out)
        & )

    else
      xy_DAtmMassDt                  = 0.0_DP
      xy_LatHeatFluxByMajCompIceSubl = 0.0_DP
    end if

    xy_DPsDt = xy_DAtmMassDt * Grav

    xy_LatHeatFluxByOtherSpc = xy_LatHeatFluxByMajCompIceSubl

    if ( FlagSublimation ) then
      ! If sublimation is considered, the melt of snow/ice is not calculated.
      xy_LatHeatFluxBySnowMelt = 0.0_DP
    else
      ! Else, the melt of snow/ice is calculated.
      call PhyImplSDHV3IceSnowPhaseChgCor(               &
        & IndexSpcH2O,                                   & ! (in)
        & xyr_Press(:,:,0),                              & ! (in)
        & xyr_HeatFlux, xy_SurfLatentHeatFlux,           & ! (in)
        & xyr_SoilHeatFlux,                              & ! (in)
        & xyr_SoilTempTransCoef,                         & ! (in)
        & xyr_RadSFlux, xyr_RadLFlux,                    & ! (in)
        & xy_DeepSubSurfHeatFlux,                        & ! (in)
        & xy_SurfTemp, xyz_SoilTemp,                     & ! (in)
        & xy_SoilMoistB, xy_SurfSnowB,                   & ! (in)
        & xy_SurfHeatCapacity,                           & ! (in)
        & xy_SoilHeatCap, xy_SoilHeatDiffCoef,           & ! (in)
        & xy_IndexCalcMethod,                            & ! (in)
        & xyra_DelRadLFlux,                              & ! (in)
        & xyz_Exner, xyr_Exner,                          & ! (in)
        & xy_SurfTempTransCoef,                          & ! (in)
        & xy_LatHeatFluxByOtherSpc,                      & ! (in)
        & xyza_TempMtx, xyz_TempVec,                     & ! (in)
        & xyaa_SurfMtx, xy_SurfRH,                       & ! (in)
        & xyaa_SoilTempMtx, xya_SoilTempVec,             & ! (in)
        & xyz_DTempDt,                                   & ! (in)
        & xy_DSurfTempDt,                                & ! (inout)
        & xyz_DSoilTempDt,                               & ! (inout)
        & xy_DSoilMoistDt, xy_DSurfSnowDt,               & ! (inout)
        & xy_LatHeatFluxBySnowMelt                       & ! (out)
        & )
    end if



    call PhyImplSDHV3SeaIceCorrection(                             &
      & xy_IndexCalcMethod,                                        & ! (in)
      & xy_SeaIceThickness,                                        & ! (in)
      & xyz_Exner, xyr_Exner,                                      & ! (in)
      & xy_SurfTemp,                                               & ! (in)
      & xyr_HeatFlux,                                              & ! (in)
      & xyr_SoilHeatFlux,                                          & ! (in)
      & xy_SurfTempTransCoef, xyr_SoilTempTransCoef,               & ! (in)
      & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux,              & ! (in)
      & xy_SurfLatentHeatFlux,                                     & ! (in)
      & xyza_TempMtx, xyz_TempVec,                                 & ! (in)
      & xyaa_SurfMtx, xy_SurfRH,                                   & ! (in)
      & xyaa_SoilTempMtx, xya_SoilTempVec,                         & ! (in)
      & xyz_DTempDt,                                               & ! (inout)
      & xy_DSurfTempDt,                                            & ! (inout)
      & xyz_DSoilTempDt,                                           & ! (inout)
      & xy_LatHeatFluxBySeaIceMelt                                 & ! (out)
      & )


    do l = -1, 1
      do k = 1, kmax
        xyza_QMixLUMtx(:,:,k,l) = xyza_QMixMtx(:,:,k,l)
      end do
    end do

    call PhyImplLUDecomp3( &
      & xyza_QMixLUMtx,    &   ! (inout)
      & imax * jmax, kmax  &   ! (in)
      & )

    do n = 1, ncmax
      do k = 1, kmax
        xyz_DelQMixLUVec(:,:,k) = xyzf_QMixVec(:,:,k,n)
      end do

      call PhyImplLUSolve3(      &
        & xyz_DelQMixLUVec,      & ! (inout)
        & xyza_QMixLUMtx,        & ! (in)
        & 1, imax * jmax , kmax  & ! (in)
        & )

      do k = 1, kmax
        xyzf_DQMixDt(:,:,k,n) = xyz_DelQMixLUVec(:,:,k) / ( 2.0_DP * DelTime )
      end do
    end do



    ! Debug routine
    !
    call PhyImplSDHV3ChkConservation(                          &
      & xy_IndexCalcMethod,                                    &
      & xy_SeaIceThickness,                                     &
      & xyr_Press, xyz_Exner, xyr_Exner,                       &
      & xy_SurfTemp,                                           &
      & xy_SurfHeatCapacity, xy_SoilHeatCap,                   &
!!$      & xyr_HeatFlux, xyrf_QMixFlux,                           &
      & xyr_HeatFlux,                                          &
      & xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux,              &
      & xyr_SoilHeatFlux,                                      &
      & xy_SurfTempTransCoef, xyr_SoilTempTransCoef,           &
      & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux,          &
      & xy_LatHeatFluxByMajCompIceSubl,                        &
      & xy_LatHeatFluxBySnowMelt, xy_LatHeatFluxBySeaIceMelt,  &
      & xy_DeepSubSurfHeatFlux,                                &
      & xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt,           &
      & xy_DSoilMoistDt, xy_DSurfSnowDt,                                      &
      & xy_DPsDt, xy_DSurfMajCompIceDt                                        &
      & )


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

  end subroutine PhyImplSDHV3TendencyHeatCore

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

  subroutine PhyImplSDHV3TendencyHeatTQCore(                   &
    & xy_IndexCalcMethod,                                      & ! (in)
    & xy_SeaIceThickness,                                      & ! (in)
    & xy_SnowFrac,                                             & ! (in)
    & xyr_HeatFlux, xyrf_QMixFlux,                             & ! (in)
    & xy_SurfSoilHeatFlux,                                     & ! (in)
    & xyr_RadSFlux, xyr_RadLFlux,                              & ! (in)
    & xy_DeepSubSurfHeatFlux,                                  & ! (in)
    & xy_SurfTemp,                                             & ! (in)
    & xy_SurfHumidCoef,                                        & ! (in)
    & xy_SurfHeatCapacity,                                     & ! (in)
    & xyra_DelRadLFlux,                                        & ! (in)
    & xyr_Press, xyz_Exner, xyr_Exner,                         & ! (in)
    & xyr_VelTransCoef, xyr_TempTransCoef,                     & ! (in)
    & xyr_QMixTransCoef,                                       & ! (in)
    & xy_SurfVelTransCoef, xy_SurfTempTransCoef,               & ! (in)
    & xy_SurfQVapTransCoef,                                    & ! (in)
    & xyz_DTempDt, xyzf_DQMixDt,                               & ! (out)
    & xy_DSurfTempDt                                           & ! (out)
    & )
    !
    ! 時間変化率の計算を行います. 
    !
    ! Calculate tendencies. 
    !

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

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

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

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

    ! 飽和比湿の算出
    ! Evaluation of saturation specific humidity
    !
    use saturate, only: &
      & xy_CalcQVapSatOnLiq, &
      & xy_CalcQVapSatOnSol, &
      & xy_CalcDQVapSatDTempOnLiq, &
      & xy_CalcDQVapSatDTempOnSol

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method
    real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              !
                              ! Sea ice thickness
    real(DP), intent(in):: xy_SnowFrac  (0:imax-1, 1:jmax)
                              !
                              ! Snow fraction
    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux

    real(DP), intent(in):: xy_SurfSoilHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面土壌熱伝導フラックス.
                              ! Soil heat conduction flux at the surface

    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux

    real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
                              ! 地中熱フラックス. 
                              ! "Deep subsurface heat flux"
                              ! Heat flux at the bottom of surface/soil layer.
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
                              ! 地表熱容量. 
                              ! Surface heat capacity

    real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(in):: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass of constituents

    real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity

    real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(out):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 質量混合比変化. 
                              ! Mass mixing ratio tendency
    real(DP), intent(out):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP):: xyza_QMixMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 質量混合比陰解行列. 
                              ! Implicit matrix about mass mixing ratio
    real(DP):: xyz_QMixVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 質量混合比陰解ベクトル. 
                              ! Implicit vector about mass mixing ratio
    real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance
    real(DP):: xy_SurfRH(0:imax-1,1:jmax)


    real(DP):: xy_SurfQVapSatOnLiq(0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface

    real(DP):: xyaa_TempQVapLUMtx (0:imax-1, 1:jmax, -kmax:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xya_DelTempQVapLUVec (0:imax-1, 1:jmax, -kmax:kmax)
                              ! $ T, Qv $ の時間変化.
                              ! Tendency of $ T $ and $ Qv $

    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               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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



!!$    if ( .not. FlagSSModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$    end if

    ! FlagBucketModel は関係ないよね?
    ! SSModel 強制にした時点で, 水蒸気は地面と分離したから. 
!!$    if ( .not. FlagBucketModel ) then
!!$      call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$    end if


    ! 飽和比湿の計算
    ! Calculate saturated specific humidity
    !
    xy_SurfQVapSatOnLiq       = &
      & xy_CalcQVapSatOnLiq      ( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSatOnSol       = &
      & xy_CalcQVapSatOnSol      ( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSat       = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq &
      & + xy_SnowFrac              * xy_SurfQVapSatOnSol
    xy_SurfDQVapSatDTempOnLiq = &
      & xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
    xy_SurfDQVapSatDTempOnSol = &
      & xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
    xy_SurfDQVapSatDTemp = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq &
      & + xy_SnowFrac              * xy_SurfDQVapSatDTempOnSol


    ! 陰解法のための行列作成
    ! Create matrices for implicit scheme
    !

    ! 鉛直拡散スキームの輸送係数から陰解行列の計算 (温度)
    ! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (temperature)
    !
    k = 1
    xyza_TempMtx(:,:,k,-1) = &
      & - CpDry * xy_SurfTempTransCoef(:,:)
    xyza_TempMtx(:,:,k, 0) =                                                          &
      & - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k  ) * xy_SurfTempTransCoef(:,:) &
      & + CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k  )
    xyza_TempMtx(:,:,k, 1) = &
      & - CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k  )

    do k = 2, kmax-1
      xyza_TempMtx(:,:,k,-1) = &
        & - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
      xyza_TempMtx(:,:,k, 0) =                                                          &
        & - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k-1)&
        & + CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k  )
      xyza_TempMtx(:,:,k, 1) = &
        & - CpDry * xyr_Exner(:,:,k  ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k  )
    end do

    k = kmax
    xyza_TempMtx(:,:,k,-1) = &
      & - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
    xyza_TempMtx(:,:,k, 0) =                                                          &
      & - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k  ) * xyr_TempTransCoef(:,:,k-1)
    xyza_TempMtx(:,:,k, 1) = 0.0_DP

    do k = 1, kmax
      xyz_TempVec(:,:,k) = - ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) )
    end do


    ! 鉛直拡散スキームの輸送係数から陰解行列の計算 (比湿)
    ! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (specific humidity)
    !

    k = 1
    xyza_QMixMtx(:,:,k,-1) =                                                  &
        & - xy_SurfHumidCoef(:,:) * xy_SurfQVapTransCoef(:,:)                 &
        &   * xy_SurfDQVapSatDTemp(:,:)
    xyza_QMixMtx(:,:,k, 0) =                                                  &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xy_SurfHumidCoef(:,:) * xy_SurfQVapTransCoef(:,:)                   &
      & + xyr_QMixTransCoef(:,:,k  )
    xyza_QMixMtx(:,:,k, 1) =                                                  &
      & - xyr_QMixTransCoef(:,:,k  )

    do k = 2, kmax-1
      xyza_QMixMtx(:,:,k,-1) =                                                  &
        & - xyr_QMixTransCoef(:,:,k-1)
      xyza_QMixMtx(:,:,k, 0) =                                                  &
        & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + xyr_QMixTransCoef(:,:,k-1)                                          &
        & + xyr_QMixTransCoef(:,:,k  )
      xyza_QMixMtx(:,:,k, 1) =                                                  &
        & - xyr_QMixTransCoef(:,:,k  )
    end do

    k = kmax
    xyza_QMixMtx(:,:,k,-1) =                                                  &
      & - xyr_QMixTransCoef(:,:,k-1)
    xyza_QMixMtx(:,:,k, 0) =                                                  &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xyr_QMixTransCoef(:,:,k-1)
    xyza_QMixMtx(:,:,k, 1) = 0.0_DP

    n = IndexH2OVap
    do k = 1, kmax
      xyz_QMixVec(:,:,k) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
    end do



    ! 地表面過程の輸送係数から陰解行列の計算
    ! Calculate implicit matrices from transfer coefficient of surface process
    !
    do i = 0, imax-1
      do j = 1, jmax
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          !    for thermal diffusion in soil
!!$          xyaa_SurfMtx(i,j,0,-1) =                                             &
!!$            &   xyr_SoilTempTransCoef(i,j,0)
          !    for water vapor diffusion in atmosphere
          xyaa_SurfMtx(i,j,0,-1) =                                             &
            & - LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)
          xyaa_SurfMtx(i,j,0, 0) =                                             &
            &   xy_SurfHeatCapacity(i,j) / ( 2.0_DP * DelTime )                    &
            & + CpDry * xy_SurfTempTransCoef(i,j)                              &
            & + xyra_DelRadLFlux(i,j,0,0)                                      &
            & + LatentHeat                                                     &
            &   * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)            &
            &   * xy_SurfDQVapSatDTemp(i,j)
          xyaa_SurfMtx(i,j,0, 1) =                                             &
            & - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1)                    &
            &   * xy_SurfTempTransCoef(i,j)                                    &
            & + xyra_DelRadLFlux(i,j,0,1)
        case ( IndexSeaIce )
          ! sea ice
!!$          xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
          !    for water vapor diffusion in atmosphere
          xyaa_SurfMtx(i,j,0,-1) =                                             &
            & - LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)
          xyaa_SurfMtx(i,j,0, 0) =                                        &
            &   SeaIceVolHeatCap * xy_SeaIceThickness(i,j)                &
            & / ( 2.0_DP * DelTime )  &
            & + CpDry * xy_SurfTempTransCoef(i,j)                         &
            & + xyra_DelRadLFlux(i,j,0,0)                                 &
            & + SeaIceThermCondCoef / xy_SeaIceThickness(i,j)             &
            & + LatentHeat                                                &
            &   * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)       &
            &   * xy_SurfDQVapSatDTemp(i,j)
          xyaa_SurfMtx(i,j,0, 1) = &
            & - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1)                    &
            &   * xy_SurfTempTransCoef(i,j)                                    &
            & + xyra_DelRadLFlux(i,j,0,1)
        case ( IndexSlabOcean )
          ! slab ocean ocean
!!$          xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
          !    for water vapor diffusion in atmosphere
          xyaa_SurfMtx(i,j,0,-1) =                                             &
            & - LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)
          xyaa_SurfMtx(i,j,0, 0) =                                             &
            &   SOHeatCapacity / ( 2.0_DP * DelTime )                              &
            & + CpDry * xy_SurfTempTransCoef(i,j)                              &
            & + xyra_DelRadLFlux(i,j,0,0)                                      &
            & + LatentHeat                                                &
            &   * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)       &
            &   * xy_SurfDQVapSatDTemp(i,j)
          xyaa_SurfMtx(i,j,0, 1) =                                             &
            & - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1)                    &
            &   * xy_SurfTempTransCoef(i,j)                                    &
            & + xyra_DelRadLFlux(i,j,0,1)
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
          xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
          xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select

      end do
    end do

    n = IndexH2OVap
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          xy_SurfRH(i,j) =                                    &
            & - xyr_RadSFlux(i,j,0)                           &
            & - xyr_RadLFlux(i,j,0)                           &
            & - xyr_HeatFlux(i,j,0)                           &
!!$            & - xy_SurfLatentHeatFlux(i,j)                    &
            & - LatentHeat * xyrf_QMixFlux(i,j,0,n)           &
!            & + xy_DeepSubSurfHeatFlux(i,j)
!!$            & + xyr_SoilHeatFlux(i,j,0)
            & + xy_SurfSoilHeatFlux(i,j)
        case ( IndexSeaIce )
          ! sea ice
          xy_SurfRH(i,j) =                                      &
            & - xyr_RadSFlux(i,j,0)                             &
            & - xyr_RadLFlux(i,j,0)                             &
            & - xyr_HeatFlux(i,j,0)                             &
!!$            & - xy_SurfLatentHeatFlux(i,j)                      &
            & - LatentHeat * xyrf_QMixFlux(i,j,0,n)             &
!              & + xy_DeepSubSurfHeatFlux(i,j)
!              & + xyr_SoilHeatFlux(i,j,0)
            & - SeaIceThermCondCoef                             &
            &   * ( xy_SurfTemp(i,j) - TempBelowSeaIce ) / xy_SeaIceThickness(i,j)
        case ( IndexSlabOcean )
          ! slab ocean
          xy_SurfRH(i,j) =                                       &
            & - xyr_RadSFlux(i,j,0)                              &
            & - xyr_RadLFlux(i,j,0)                              &
            & - xyr_HeatFlux(i,j,0)                              &
!!$            & - xy_SurfLatentHeatFlux(i,j)                       !&
            & - LatentHeat * xyrf_QMixFlux(i,j,0,n)              !&
!              & + xy_DeepSubSurfHeatFlux(i,j)
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xy_SurfRH(i,j) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do


    ! 温度と比湿の計算
    ! Calculate temperature and specific humidity
    !

    do l = -1, 1
      do k = 1, kmax
        xyaa_TempQVapLUMtx(:,:,-k,-l) = xyza_QMixMtx(:,:,k,l)
      end do
      k = 0
      xyaa_TempQVapLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
      do k = 1, kmax
        xyaa_TempQVapLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
      end do
    end do

    call PhyImplLUDecomp3(              &
      & xyaa_TempQVapLUMtx,             & ! (inout)
      & imax * jmax, kmax + 1 + kmax    & ! (in)
      )

    do k = 1, kmax
      xya_DelTempQVapLUVec(:,:,-k) = xyz_QMixVec(:,:,k)
    end do
    k = 0
    xya_DelTempQVapLUVec(:,:,k) = xy_SurfRH
    do k = 1, kmax
      xya_DelTempQVapLUVec(:,:,k) = xyz_TempVec(:,:,k)
    end do

    call PhyImplLUSolve3(                   &
      & xya_DelTempQVapLUVec,               & ! (inout)
      & xyaa_TempQVapLUMtx,                 & ! (in)
      & 1, imax * jmax , kmax + 1 + kmax    & ! (in)
      & )

    n = IndexH2OVap
    do k = 1, kmax
      xyzf_DQMixDt(:,:,k,n) = xya_DelTempQVapLUVec(:,:,-k) / ( 2.0_DP * DelTime )
    end do
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          xy_DSurfTempDt(i,j) = xya_DelTempQVapLUVec(i,j,0) / ( 2.0_DP * DelTime )
        case ( IndexSeaIce )
          ! sea ice
          xy_DSurfTempDt(i,j) = xya_DelTempQVapLUVec(i,j,0) / ( 2.0_DP * DelTime )
        case ( IndexSlabOcean )
          ! slab ocean
          xy_DSurfTempDt(i,j) = xya_DelTempQVapLUVec(i,j,0) / ( 2.0_DP * DelTime )
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xy_DSurfTempDt(i,j) = 0.
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    do k = 1, kmax
      xyz_DTempDt(:,:,k) = xya_DelTempQVapLUVec(:,:,k) / ( 2.0_DP * DelTime )
    end do


    call PhyImplSDHV3ChkConservationTQ(                   &
      & xy_IndexCalcMethod,                                     &
      & xy_SeaIceThickness,                                     &
      & xy_SnowFrac,                                            &
      & xyr_Press, xyz_Exner, xyr_Exner,                        &
      & xy_SurfTemp,                                            &
      & xyr_HeatFlux, xyrf_QMixFlux,                            &
      & xy_SurfSoilHeatFlux,                                    &
      & xy_SurfTempTransCoef, xy_SurfQVapTransCoef,             &
      & xy_SurfHumidCoef, xy_SurfHeatCapacity,                  &
      & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux,           &
      & xy_DeepSubSurfHeatFlux,                                 &
      & xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt               &
      & )


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

  end subroutine PhyImplSDHV3TendencyHeatTQCore

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

  subroutine PhyImplSDHV3ChkConservation(                     &
    & xy_IndexCalcMethod,                                     &
    & xy_SeaIceThickness,                                     &
    & xyr_Press, xyz_Exner, xyr_Exner,                        &
    & xy_SurfTemp,                                            &
    & xy_SurfHeatCapacity, xy_SoilHeatCap,                    &
!!$    & xyr_HeatFlux, xyrf_QMixFlux,                            &
    & xyr_HeatFlux,                                           &
    & xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux,               &
    & xyr_SoilHeatFlux,                                       &
    & xy_SurfTempTransCoef, xyr_SoilTempTransCoef,            &
    & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux,           &
    & xy_LatHeatFluxByMajCompIceSubl,                         &
    & xy_LatHeatFluxBySnowMelt, xy_LatHeatFluxBySeaIceMelt,   &
    & xy_DeepSubSurfHeatFlux,                                 &
    & xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt,            &
    & xy_DSoilMoistDt, xy_DSurfSnowDt,                                       &
    & xy_DPsDt, xy_DSurfMajCompIceDt                                         &
    & )
    !
    ! 
    !
    ! A part of conservation of energy is checked.
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_SSDepth,   &    ! subsurface grid on interface of layer
      & z_SSDepth         ! subsurface grid at midpoint of layer

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

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

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

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method

    real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              !
                              ! Sea ice thickness

    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature

    real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
                              ! 地表熱容量. 
                              ! Surface heat capacity
    real(DP), intent(in):: xy_SoilHeatCap (0:imax-1, 1:jmax)
                              ! 土壌熱容量 (J K-1 kg-1)
                              ! Specific heat of soil (J K-1 kg-1)

    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
!!$    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
!!$                              ! 比湿フラックス. 
!!$                              ! Specific humidity flux
    real(DP), intent(in):: xy_SurfH2OVapFlux(0:imax-1, 1:jmax)
                              ! 惑星表面水蒸気フラックス.
                              ! Water vapor flux at the surface
    real(DP), intent(in):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface
    real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
                              ! 土壌中の熱フラックス (W m-2)
                              ! Heat flux in sub-surface soil (W m-2)
    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
                              ! 輸送係数：土壌温度.
                              ! Transfer coefficient: soil temperature

    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in ):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by major component ice sublimation
                              ! (variable only for debug)
    real(DP), intent(in ):: xy_LatHeatFluxBySnowMelt      (0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by melt (variable only for debug)
    real(DP), intent(in ):: xy_LatHeatFluxBySeaIceMelt    (0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by sea ice melt (variable only for debug)
    real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
                              ! 地中熱フラックス. 
                              ! "Deep subsurface heat flux"
                              ! Heat flux at the bottom of surface/soil layer.

    real(DP), intent(in ):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(in ):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 質量混合比変化. 
                              ! Mass mixing ratio tendency
    real(DP), intent(in ):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)
    real(DP), intent(in ):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
                              ! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
                              ! Temperature tendency (K s-1)
    real(DP), intent(in ):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
                              ! 土壌温度時間変化率 (kg m-2 s-1)
                              ! Soil temperature tendency (kg m-2 s-1)
    real(DP), intent(in ):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
                              ! 積雪率時間変化率 (kg m-2 s-1)
                              ! Surface snow amount tendency (kg m-2 s-1)

    real(DP), intent(in ):: xy_DPsDt            (0:imax-1, 1:jmax)
    real(DP), intent(in ):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)


    ! 作業変数
    ! Work variables
    !

    real(DP) :: xy_SurfRadSFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfRadLFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSensHeatFlux(0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
    real(DP) :: xy_SeaIceHeatCondFlux  (0:imax-1, 1:jmax)

    real(DP) :: xy_Residual            (0:imax-1, 1:jmax)
    real(DP) :: xy_SumAtmRate          (0:imax-1, 1:jmax)

    real(DP) :: MaxResidual

    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
    !

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


    xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
    xy_SurfRadLFlux = xyr_RadLFlux(:,:,0)                                           &
      &   + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime )        &
      &   + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
    xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0)                                &
      & - xyr_SoilTempTransCoef(:,:,0)                                               &
      &   * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0_DP * DelTime )
    xy_SurfSensHeatFlux =                                                       &
      & xyr_HeatFlux(:,:,0)                                                     &
      &   - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef                     &
      &     * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1)                           &
      &       - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
    xy_SeaIceHeatCondFlux =                                                           &
      & - SeaIceThermCondCoef                                                         &
      &   * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) &
      &   / xy_SeaIceThickness


    !-----
    ! Atmospheric heating
    !
    xy_SumAtmRate = 0.0_DP
    do k = kmax, 1, -1
      xy_SumAtmRate = xy_SumAtmRate                                  &
        & + CpDry * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav &
        &     * xyz_DTempDt(:,:,k)
    end do
    !
    xy_Residual =                          &
      & - ( 0.0_DP - xy_SurfSensHeatFlux ) &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )

!!$        select case ( xy_IndexCalcMethod(i,j) )
!!$        case ( IndexSeaIce )
!!$        case default
!!$          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$        end select

      end do
    end do
    if ( MaxResidual > 1.0d-10 ) then
      call MessageNotify( 'M', module_name, &
        & 'Atm. sensible heating res. : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Land surface
    !
    xy_SumAtmRate = xy_SurfHeatCapacity * xy_DSurfTempDt
    !
    xy_Residual =                        &
      &   xy_SurfRadSFlux                &
      & + xy_SurfRadLFlux                &
      & + xy_SurfSensHeatFlux            &
      & + xy_SurfLatentHeatFlux          &
      & - xy_SurfSoilHeatCondFlux        &
      & + xy_LatHeatFluxByMajCompIceSubl &
      & + xy_LatHeatFluxBySnowMelt       &
      & + xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexSeaIce )
          ! sea ice
        case ( IndexSlabOcean )
          ! slab ocean
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Land surf. heat budget res.: %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Soil heating
    !
    xy_SumAtmRate = 0.0_DP
    do k = 1, kslmax
      xy_SumAtmRate = xy_SumAtmRate                            &
        & + xy_SoilHeatCap * ( r_SSDepth(k-1) - r_SSDepth(k) ) &
        &     * xyz_DSoilTempDt(:,:,k)
    end do
    !
    xy_Residual =                                               &
      & - ( xy_SurfSoilHeatCondFlux - xy_DeepSubSurfHeatFlux )  &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexSeaIce )
          ! sea ice
        case ( IndexSlabOcean )
          ! slab ocean
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Soil heating res.          : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Slab ocean heating
    !
    xy_SumAtmRate = SOHeatCapacity * xy_DSurfTempDt
    !
    xy_Residual =                       &
      & - (                             &
      &       xy_SurfRadSFlux           &
      &     + xy_SurfRadLFlux           &
      &     + xy_SurfSensHeatFlux       &
      &     + xy_SurfLatentHeatFlux     &
      &   )                             &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
        case ( IndexSeaIce )
          ! sea ice
        case ( IndexSlabOcean )
          ! slab ocean
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Slab ocean heating res.    : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Sea ice heating
    !
    xy_SumAtmRate = SeaIceVolHeatCap * xy_SeaIceThickness * xy_DSurfTempDt
    !
    xy_Residual =                         &
      & - (                               &
      &     (                             &
      &         xy_SurfRadSFlux           &
      &       + xy_SurfRadLFlux           &
      &       + xy_SurfSensHeatFlux       &
      &       + xy_SurfLatentHeatFlux     &
      &     )                             &
      &     - xy_SeaIceHeatCondFlux       &
      &   )                               &
      & - xy_LatHeatFluxBySeaIceMelt      &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
        case ( IndexSeaIce )
          ! sea ice
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexSlabOcean )
          ! slab ocean
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Sea ice heating res.       : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Atmospheric moistening
    !
    xy_SumAtmRate = 0.0_DP
    do k = kmax, 1, -1
      xy_SumAtmRate = xy_SumAtmRate                          &
        & + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav &
        &     * xyzf_DQMixDt(:,:,k,IndexH2OVap)
    end do
    !
    xy_Residual =                                       &
!!$      & - ( 0.0_DP - xyrf_QMixFlux(:,:,0,IndexH2OVap) ) &
      & - ( 0.0_DP - xy_SurfH2OVapFlux ) &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Atm. moistening res.       : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Land water budget
    !
    xy_SumAtmRate = 0.0_DP
    do k = kmax, 1, -1
      xy_SumAtmRate = xy_SumAtmRate                          &
        & + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav &
        &     * xyzf_DQMixDt(:,:,k,IndexH2OVap)
    end do
    !
    xy_Residual = xy_DSoilMoistDt + xy_DSurfSnowDt + xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexSeaIce )
          ! sea ice
        case ( IndexSlabOcean )
          ! slab ocean
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Land water budget res.     : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Atmospheric mass budget
    !
    xy_Residual = xy_DPsDt / Grav + xy_DSurfMajCompIceDt
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Atm. mass budget res.      : %f.', d = (/ MaxResidual /) )
    end if


  end subroutine PhyImplSDHV3ChkConservation

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

  subroutine PhyImplSDHV3ChkConservationTQ(                   &
    & xy_IndexCalcMethod,                                     &
    & xy_SeaIceThickness,                                     &
    & xy_SnowFrac,                                            &
    & xyr_Press, xyz_Exner, xyr_Exner,                        &
    & xy_SurfTemp,                                            &
    & xyr_HeatFlux, xyrf_QMixFlux,                            &
    & xy_SurfSoilHeatFlux,                                    &
    & xy_SurfTempTransCoef, xy_SurfQVapTransCoef,             &
    & xy_SurfHumidCoef, xy_SurfHeatCapacity,                  &
    & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux,           &
    & xy_DeepSubSurfHeatFlux,                                 &
    & xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt               &
    & )
    !
    ! 
    !
    ! A part of conservation of energy is checked.
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_SSDepth,   &    ! subsurface grid on interface of layer
      & z_SSDepth         ! subsurface grid at midpoint of layer

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

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

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

    ! 飽和比湿の算出
    ! Evaluation of saturation specific humidity
    !
    use saturate, only: &
      & xy_CalcQVapSatOnLiq, &
      & xy_CalcQVapSatOnSol, &
      & xy_CalcDQVapSatDTempOnLiq, &
      & xy_CalcDQVapSatDTempOnSol

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method
    real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              !
                              ! Sea ice thickness
    real(DP), intent(in):: xy_SnowFrac(0:imax-1, 1:jmax)
                              !
                              ! Snow fraction
    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature

    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(in):: xy_SurfSoilHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面土壌熱伝導フラックス.
                              ! Soil heat conduction flux at the surface
    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity
    real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
                              ! 地表熱容量. 
                              ! Surface heat capacity

    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
                              ! 地中熱フラックス. 
                              ! "Deep subsurface heat flux"
                              ! Heat flux at the bottom of surface/soil layer.

    real(DP), intent(in ):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(in ):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 質量混合比変化. 
                              ! Mass mixing ratio tendency
    real(DP), intent(in ):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)


    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_SurfQVapSatOnLiq(0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface

    real(DP) :: xy_SurfRadSFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfRadLFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSensHeatFlux    (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfH2OVapFlux      (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfLatentHeatFlux  (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
    real(DP) :: xy_SeaIceHeatCondFlux  (0:imax-1, 1:jmax)

    real(DP) :: xy_Residual            (0:imax-1, 1:jmax)
    real(DP) :: xy_SumAtmRate          (0:imax-1, 1:jmax)

    real(DP) :: MaxResidual

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


    ! 実行文 ; Executable statement
    !

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


    ! 飽和比湿の計算
    ! Calculate saturated specific humidity
    !
    xy_SurfQVapSatOnLiq       = &
      & xy_CalcQVapSatOnLiq      ( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSatOnSol       = &
      & xy_CalcQVapSatOnSol      ( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSat       = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq &
      & + xy_SnowFrac              * xy_SurfQVapSatOnSol
    xy_SurfDQVapSatDTempOnLiq = &
      & xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
    xy_SurfDQVapSatDTempOnSol = &
      & xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
    xy_SurfDQVapSatDTemp = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq &
      & + xy_SnowFrac              * xy_SurfDQVapSatDTempOnSol


    xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
    xy_SurfRadLFlux = xyr_RadLFlux(:,:,0)                                           &
      &   + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime )        &
      &   + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
    xy_SurfSoilHeatCondFlux = xy_SurfSoilHeatFlux
    xy_SurfSensHeatFlux =                                                       &
      & xyr_HeatFlux(:,:,0)                                                     &
      &   - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef                     &
      &     * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1)                           &
      &       - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
    xy_SeaIceHeatCondFlux =                                                           &
      & - SeaIceThermCondCoef                                                         &
      &   * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) &
      &   / xy_SeaIceThickness

    n = IndexH2OVap
    xy_SurfH2OVapFlux = xyrf_QMixFlux(:,:,0,n)               &
      & - xy_SurfHumidCoef * xy_SurfQVapTransCoef            &
      &   * ( xyzf_DQMixDt(:,:,1,n)                          &
      &     - xy_SurfDQVapSatDTemp * xy_DSurfTempDt )        &
      &   * 2.0_DP * DelTime
    xy_SurfLatentHeatFlux = LatentHeat * xy_SurfH2OVapFlux


    !-----
    ! Atmospheric heating
    !
    xy_SumAtmRate = 0.0_DP
    do k = kmax, 1, -1
      xy_SumAtmRate = xy_SumAtmRate                                  &
        & + CpDry * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav &
        &     * xyz_DTempDt(:,:,k)
    end do
    !
    xy_Residual =                          &
      & - ( 0.0_DP - xy_SurfSensHeatFlux ) &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )

!!$        select case ( xy_IndexCalcMethod(i,j) )
!!$        case ( IndexSeaIce )
!!$        case default
!!$          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$        end select

      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Atm. sensible heating res. : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Land surface
    !
    xy_SumAtmRate = xy_SurfHeatCapacity * xy_DSurfTempDt
    !
    xy_Residual =                        &
      &   xy_SurfRadSFlux                &
      & + xy_SurfRadLFlux                &
      & + xy_SurfSensHeatFlux            &
      & + xy_SurfLatentHeatFlux          &
      & - xy_SurfSoilHeatCondFlux        &
      & + xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexSeaIce )
          ! sea ice
        case ( IndexSlabOcean )
          ! slab ocean
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Land surf. heat budget res.: %f.', d = (/ MaxResidual /) )
    end if


!!$    !-----
!!$    ! Soil heating
!!$    !
!!$    xy_SumAtmRate = 0.0_DP
!!$    do k = 1, kslmax
!!$      xy_SumAtmRate = xy_SumAtmRate                            &
!!$        & + xy_SoilHeatCap * ( r_SSDepth(k-1) - r_SSDepth(k) ) &
!!$        &     * xyz_DSoilTempDt(:,:,k)
!!$    end do
!!$    !
!!$    xy_Residual =                                               &
!!$      & - ( xy_SurfSoilHeatCondFlux - xy_DeepSubSurfHeatFlux )  &
!!$      & - xy_SumAtmRate
!!$    !
!!$    MaxResidual = 0.0_DP
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        select case ( xy_IndexCalcMethod(i,j) )
!!$        case ( IndexLand )
!!$          ! land
!!$          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$        case ( IndexSeaIce )
!!$          ! sea ice
!!$        case ( IndexSlabOcean )
!!$          ! slab ocean
!!$        case ( IndexOceanPresSST )
!!$          ! open ocean
!!$        case default
!!$          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
!!$        end select
!!$      end do
!!$    end do
!!$    if ( MaxResidual > 1.0d-10 ) then
!!$      call MessageNotify( 'M', module_name, &
!!$        & 'Soil heating res.          : %f.', d = (/ MaxResidual /) )
!!$    end if


    !-----
    ! Slab ocean heating
    !
    xy_SumAtmRate = SOHeatCapacity * xy_DSurfTempDt
    !
    xy_Residual =                       &
      & - (                             &
      &       xy_SurfRadSFlux           &
      &     + xy_SurfRadLFlux           &
      &     + xy_SurfSensHeatFlux       &
      &     + xy_SurfLatentHeatFlux     &
      &   )                             &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
        case ( IndexSeaIce )
          ! sea ice
        case ( IndexSlabOcean )
          ! slab ocean
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Slab ocean heating res.    : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Sea ice heating
    !
    xy_SumAtmRate = SeaIceVolHeatCap * xy_SeaIceThickness * xy_DSurfTempDt
    !
    xy_Residual =                         &
      & - (                               &
      &     (                             &
      &         xy_SurfRadSFlux           &
      &       + xy_SurfRadLFlux           &
      &       + xy_SurfSensHeatFlux       &
      &       + xy_SurfLatentHeatFlux     &
      &     )                             &
      &     - xy_SeaIceHeatCondFlux       &
      &   )                               &
!!$      & - xy_LatHeatFluxBySeaIceMelt      &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
        case ( IndexSeaIce )
          ! sea ice
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexSlabOcean )
          ! slab ocean
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Sea ice heating res.       : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Atmospheric moistening
    !
    xy_SumAtmRate = 0.0_DP
    do k = kmax, 1, -1
      xy_SumAtmRate = xy_SumAtmRate                          &
        & + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav &
        &     * xyzf_DQMixDt(:,:,k,IndexH2OVap)
    end do
    !
    xy_Residual =                                       &
!!$      & - ( 0.0_DP - xyrf_QMixFlux(:,:,0,IndexH2OVap) ) &
      & - ( 0.0_DP - xy_SurfH2OVapFlux ) &
      & - xy_SumAtmRate
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Atm. moistening res.       : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Land water budget
    !
    xy_SumAtmRate = 0.0_DP
    do k = kmax, 1, -1
      xy_SumAtmRate = xy_SumAtmRate                          &
        & + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav &
        &     * xyzf_DQMixDt(:,:,k,IndexH2OVap)
    end do
    !
!!$    xy_Residual = xy_DSoilMoistDt + xy_DSurfSnowDt + xy_SumAtmRate
    xy_Residual = xy_SumAtmRate - xy_SurfH2OVapFlux
    !
    MaxResidual = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
        case ( IndexSeaIce )
          ! sea ice
        case ( IndexSlabOcean )
          ! slab ocean
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    if ( MaxResidual > 1.0e-10_DP ) then
      call MessageNotify( 'M', module_name, &
        & 'Land water budget res.     : %f.', d = (/ MaxResidual /) )
    end if


!!$    !-----
!!$    ! Atmospheric mass budget
!!$    !
!!$    xy_Residual = xy_DPsDt / Grav + xy_DSurfMajCompIceDt
!!$    !
!!$    MaxResidual = 0.0_DP
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$      end do
!!$    end do
!!$    if ( MaxResidual > 1.0d-10 ) then
!!$      call MessageNotify( 'M', module_name, &
!!$        & 'Atm. mass budget res.      : %f.', d = (/ MaxResidual /) )
!!$    end if


  end subroutine PhyImplSDHV3ChkConservationTQ

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

!!$  subroutine PhyImplSDHV3IceSnowPhaseChgCor(       &
  subroutine OLD_PhyImplSDHV3IceSnowPhase(       &
    & IndexSpc,                                      & ! (in)
    & xy_Ps,                                         & ! (in)
    & xyr_HeatFlux, xy_SurfLatentHeatFlux,           & ! (in)
    & xyr_SoilHeatFlux,                              & ! (in)
    & xyr_SoilTempTransCoef,                         & ! (in)
    & xyr_RadSFlux, xyr_RadLFlux,                    & ! (in)
    & xy_DeepSubSurfHeatFlux,                        & ! (in)
    & xy_SurfTemp, xyz_SoilTemp,                     & ! (in)
    & xy_SurfLiqB, xy_SurfSolB,                      & ! (in)
    & xy_SurfHeatCapacity,                           & ! (in)
    & xy_SoilHeatCap, xy_SoilHeatDiffCoef,           & ! (in)
    & xy_IndexCalcMethod,                            & ! (in)
    & xyra_DelRadLFlux,                              & ! (in)
    & xyz_Exner, xyr_Exner,                          & ! (in)
    & xy_SurfTempTransCoef,                          & ! (in)
    & xy_LatHeatFluxByOtherSpc,                      & ! (in)
    & xyza_ArgTempMtx, xyz_ArgTempVec,               & ! (in)
    & xyaa_ArgSurfMtx, xy_ArgSurfRH,                 & ! (in)
    & xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec,       & ! (in)
    & xyz_DTempDt,                                   & ! (in)
    & xy_DSurfTempDt,                                & ! (inout)
    & xyz_DSoilTempDt,                               & ! (inout)
    & xy_DSurfLiqDt, xy_DSurfSolDt,                  & ! (inout)
    & xy_LatHeatFluxBySnowMelt                       & ! (out)
    & )
    !
    ! 融雪による時間変化率の修正を行います. 
    !
    ! Correction of tendencies due to melt of snow. 
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_SSDepth,   &    ! subsurface grid on interface of layer
      & z_SSDepth         ! subsurface grid at midpoint of layer

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & CpDry,      &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & LatentHeatFusion
                              ! $ L $ [J kg-1] . 
                              ! 融解の潜熱. 
                              ! Latent heat of fusion

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

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


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

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

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


    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: IndexSpc
    real(DP), intent(in):: xy_Ps(0:imax-1, 1:jmax)
                              ! 
                              ! Surface pressure

    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux

    real(DP), intent(in):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface

    real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
                              ! 土壌中の熱フラックス (W m-2)
                              ! Heat flux in sub-surface soil (W m-2)
    real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
                              ! 輸送係数：土壌温度.
                              ! Transfer coefficient: soil temperature


    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux

    real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
                              ! 地中熱フラックス. 
                              ! "Deep subsurface heat flux"
                              ! Heat flux at the bottom of surface/soil layer.
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xyz_SoilTemp (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度 (K)
                              ! Soil temperature (K)

    real(DP), intent(in):: xy_SurfLiqB (0:imax-1, 1:jmax)
                              ! 
                              ! Surface liquid amount
    real(DP), intent(in):: xy_SurfSolB (0:imax-1, 1:jmax)
                              ! 積雪量.
                              ! Surface snow amount.

    real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
                              ! 地表熱容量. 
                              ! Surface heat capacity
    real(DP), intent(in ):: xy_SoilHeatCap (0:imax-1, 1:jmax)
                              ! 土壌熱容量 (J K-1 kg-1)
                              ! Specific heat of soil (J K-1 kg-1)
    real(DP), intent(in ):: xy_SoilHeatDiffCoef (0:imax-1, 1:jmax)
                              ! 土壌熱伝導係数 (J m-3 K-1)
                              ! Heat conduction coefficient of soil (J m-3 K-1)

    integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method

    real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature

    real(DP), intent(in   ):: xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP), intent(in   ):: xyz_ArgTempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP), intent(in   ):: xyaa_ArgSurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance
    real(DP), intent(in   ):: xy_ArgSurfRH(0:imax-1,1:jmax)


    real(DP), intent(in   ):: xyaa_ArgSoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP), intent(in   ):: xya_ArgSoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature

    real(DP), intent(in   ):: xy_LatHeatFluxByOtherSpc(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by other specie

    real(DP), intent(inout):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(inout):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)
    real(DP), intent(inout):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
                              ! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
                              ! Temperature tendency (K s-1)

    real(DP), intent(inout):: xy_DSurfLiqDt (0:imax-1, 1:jmax)
                              ! 土壌温度時間変化率 (kg m-2 s-1)
                              ! Soil temperature tendency (kg m-2 s-1)
    real(DP), intent(inout):: xy_DSurfSolDt (0:imax-1, 1:jmax)
                              ! 積雪率時間変化率 (kg m-2 s-1)
                              ! Surface snow amount tendency (kg m-2 s-1)
    real(DP), intent(out  ):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by melt (variable only for debug)

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

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

    logical :: xy_FlagCalc(0:imax-1, 1:jmax)
    integer :: xy_IndexMeltOrFreeze(0:imax-1, 1:jmax)
    integer, parameter :: IndexOthers = 0
    integer, parameter :: IndexMelt   = 1
    integer, parameter :: IndexFreeze = 2

    real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance
    real(DP):: xy_SurfRH(0:imax-1,1:jmax)


    real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature

    real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
                              ! $ T, Tg $ の時間変化.
                              ! Tendency of $ T $ and $ Tg |

    real(DP):: LatentHeatLocal
    real(DP):: LatentHeatFluxByMelt
    real(DP):: SenHeatFluxA
    real(DP):: LatHeatFluxA
    real(DP):: CondHeatFluxA
    real(DP):: ValueAlpha
    real(DP):: SurfTempATentative
    real(DP):: SoilTempATentative
    real(DP):: SurfLiqATentative
    real(DP):: xy_SurfLiqATentativeSave(0:imax-1, 1:jmax)
    real(DP):: SurfSolATentative
    real(DP):: xy_SurfSolATentativeSave(0:imax-1, 1:jmax)
    real(DP):: DelSurfSol

    real(DP) :: xy_TempMajCompCond(0:imax-1, 1:jmax)
    real(DP) :: SurfMajCompIceATentative

    real(DP) :: xy_SurfRadSFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfRadLFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSensHeatFlux    (0:imax-1, 1:jmax)
    real(DP) :: xy_SeaIceHeatCondFlux  (0:imax-1, 1:jmax)
    real(DP) :: xy_HeatingTendency     (0:imax-1, 1:jmax)

    real(DP) :: LatHeatFluxBySnowMelt

    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


    ! 実行文 ; Executable statement
    !

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

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


    !
    ! check flag of snow melt
    !
    if ( ( .not. FlagBucketModel ) .or. ( .not. FlagSnow ) ) then
      xy_LatHeatFluxBySnowMelt = 0.0_DP
      return
    end if


    if ( kslmax == 0 ) then

      call MessageNotify( 'E', module_name, &
        & 'kslmax <= 0 in PhyImplSDHV3IceSnowPhaseChgCor.' )

    else

      xy_DSurfLiqDtSave = xy_DSurfLiqDt
      xy_DSurfSolDtSave = xy_DSurfSolDt


      select case ( IndexSpc )
      case ( IndexSpcMajComp )
        call SaturateMajorCompCondTemp( &
          & xy_Ps,                      & ! (in)
          & xy_TempCond                 & ! (out)
          & )
        LatentHeatLocal = SaturateMajorCompInqLatentHeat()
      case ( IndexSpcH2O )
        xy_TempCond     = TempCondWater
        LatentHeatLocal = LatentHeatFusion
      case default
        call MessageNotify( 'E', module_name, &
          & 'Undefined IndexSpc, %d.', i = (/ IndexSpc /) )
      end select


      xy_SurfLiqATentativeSave = xy_SurfLiqB + xy_DSurfLiqDt * ( 2.0_DP * DelTime )
      xy_SurfSolATentativeSave = xy_SurfSolB + xy_DSurfSolDt * ( 2.0_DP * DelTime )

      !----------
      ! A case that a part of snow/ice melt or soil moisture freeze
      !----------

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

          select case ( xy_IndexCalcMethod(i,j) )
          case ( IndexLand )

            SurfTempATentative = xy_SurfTemp(i,j)          &
              & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
            SurfSolATentative = xy_SurfSolATentativeSave(i,j)

            if ( &
              & ( SurfSolATentative  > 0.0_DP           ) .and. &
              & ( SurfTempATentative > xy_TempCond(i,j) )       &
              & ) then
              xy_FlagCalc         (i,j) = .true.
            else
              xy_FlagCalc         (i,j) = .false.
            end if

          case default
            xy_FlagCalc         (i,j) = .false.
          end select

        end do
      end do

      ! Freeze
      select case ( IndexSpc )
      case ( IndexSpcMajComp )

        do j = 1, jmax
          do i = 0, imax-1
            select case ( xy_IndexCalcMethod(i,j) )
            case ( IndexLand )
              SurfTempATentative = xy_SurfTemp(i,j)          &
                & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
              if ( SurfTempATentative < xy_TempCond(i,j) ) then
                xy_FlagCalc         (i,j) = .true.
              end if
            end select
          end do
        end do

!!$      case ( IndexSpcH2O )
!!$
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            select case ( xy_IndexCalcMethod(i,j) )
!!$            case ( IndexLand )
!!$              SurfTempATentative = xy_SurfTemp(i,j)          &
!!$                & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
!!$              SurfLiqATentative = xy_SurfLiqATentativeSave(i,j)
!!$              if ( &
!!$                & ( SurfLiqATentative  > 0.0_DP           ) .and. &
!!$                & ( SurfTempATentative < xy_TempCond(i,j) )       &
!!$                & ) then
!!$                xy_FlagCalc         (i,j) = .true.
!!$              end if
!!$            end select
!!$
!!$          end do
!!$        end do
!!$
      end select


      xyza_TempMtx     = xyza_ArgTempMtx
      xyz_TempVec = xyz_ArgTempVec
      !
      xyaa_SurfMtx     = xyaa_ArgSurfMtx
      xy_SurfRH        = xy_ArgSurfRH
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
            xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
            xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
            xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
            xy_SurfRH   (i,j)      = xy_TempCond(i,j) - xy_SurfTemp(i,j)
          end if
        end do
      end do
      !
      xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
      xya_SoilTempVec = xya_ArgSoilTempVec


      ! 温度の計算
      ! Calculate temperature and specific humidity
      !
      do l = -1, 1
        do k = 1, kslmax
          xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
        end do
        k = 0
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
        do k = 1, kmax
          xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
        end do
      end do

      call PhyImplLUDecomp3(              &
        & xyaa_TempSoilTempLUMtx,         & ! (inout)
        & imax * jmax, kmax + 1 + kslmax  & ! (in)
        )

      do k = 1, kslmax
        xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
      end do
      k = 0
      xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
      do k = 1, kmax
        xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
      end do

      call PhyImplLUSolve3(                   &
        & xya_DelTempSoilTempLUVec,           & ! (inout)
        & xyaa_TempSoilTempLUMtx,             & ! (in)
        & 1, imax * jmax , kmax + 1 + kslmax  & ! (in)
        & )


      do k = 1, kslmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              select case ( xy_IndexCalcMethod(i,j) )
              case ( IndexLand )
                xyz_DSoilTempDt(i,j,k) = &
                  & xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
              case default
                xyz_DSoilTempDt(i,j,k) = 0.0_DP
              end select
            end if
          end do
        end do
      end do
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
            select case ( xy_IndexCalcMethod(i,j) )
            case ( IndexLand )
              ! land
              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSeaIce )
!!$              ! sea ice
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSlabOcean )
!!$              ! slab ocean
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexOceanPresSST )
!!$              ! open ocean
!!$              xy_DSurfTempDt(i,j) = 0.0_DP
            case default
              call MessageNotify( 'E', module_name, 'Unexpected Error.' )
            end select
          end if
        end do
      end do
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
            end if
          end do
        end do
      end do


      !----------
      ! Surface fluxes used below
      !----------
      xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
      xy_SurfRadLFlux = xyr_RadLFlux(:,:,0)                                           &
        &   + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime )        &
        &   + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
      xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0)                                &
        & - xyr_SoilTempTransCoef(:,:,0)                                               &
        &   * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0d0 * DelTime )
      xy_SurfSensHeatFlux =                                                       &
        & xyr_HeatFlux(:,:,0)                                                     &
        &   - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef                     &
        &     * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1)                           &
        &       - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )

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

          if ( xy_FlagCalc(i,j) ) then

            xy_LatHeatFluxBySnowMelt(i,j) =           &
              & - xy_SurfRadSFlux(i,j)                &
              & - xy_SurfRadLFlux(i,j)                &
              & - xy_SurfSensHeatFlux(i,j)            &
              & - xy_SurfLatentHeatFlux(i,j)          &
              & + xy_SurfSoilHeatCondFlux(i,j)        &
              & - xy_LatHeatFluxByOtherSpc(i,j)       &
              & - xy_SurfHeatCapacity(i,j) * xy_DSurfTempDt(i,j)

            xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
              & + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal

!!$            if ( xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * ( 2.0_DP * DelTime ) < 0.0_DP ) then
!!$              call MessageNotify( 'M', module_name, &
!!$                & 'Surface snow amount is negative %f, %f.', &
!!$                & d = (/ xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * ( 2.0_DP * DelTime ), xy_SurfSnowB(i,j) /) )
!!$            end if

          else
            xy_LatHeatFluxBySnowMelt(i,j) = 0.0_DP
          end if

        end do
      end do

      !----------
      ! A case that all snow melt or soil moisture freeze
      !----------

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

          select case ( xy_IndexCalcMethod(i,j) )
          case ( IndexLand )

            if ( xy_FlagCalc(i,j) ) then
              SurfSolATentative = xy_SurfSolB(i,j) &
                & + xy_DSurfSolDt(i,j) * 2.0_DP * DelTime
              if ( SurfSolATentative < 0.0_DP ) then
                xy_FlagCalc(i,j) = .true.
                xy_IndexMeltOrFreeze(i,j) = IndexMelt
              else
                xy_FlagCalc(i,j) = .false.
                xy_IndexMeltOrFreeze(i,j) = IndexOthers
              end if

!!$              select case ( IndexSpc )
!!$              case ( IndexSpcH2O )
!!$                SurfLiqATentative = xy_SurfLiqB(i,j) &
!!$                  & + xy_DSurfLiqDt(i,j) * 2.0_DP * DelTime
!!$                if ( SurfLiqATentative < 0.0_DP ) then
!!$                  xy_FlagCalc(i,j) = .true.
!!$                  xy_IndexMeltOrFreeze(i,j) = IndexFreeze
!!$                end if
!!$              end select

            else
              xy_FlagCalc(i,j) = .false.
              xy_IndexMeltOrFreeze(i,j) = IndexOthers
            end if

          case default
            xy_FlagCalc(i,j) = .false.
            xy_IndexMeltOrFreeze(i,j) = IndexOthers
          end select

        end do
      end do


      xyza_TempMtx     = xyza_ArgTempMtx
      xyz_TempVec      = xyz_ArgTempVec
      !
      xyaa_SurfMtx     = xyaa_ArgSurfMtx
      xy_SurfRH        = xy_ArgSurfRH
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
            select case ( xy_IndexMeltOrFreeze(i,j) )
            case ( IndexMelt )
              ! all ice/snow melt
              DelSurfSol =   xy_SurfSolATentativeSave(i,j)
            case ( IndexFreeze )
              ! all soil moisture freeze (= negative melt of ice/snow)
              DelSurfSol = - xy_SurfLiqATentativeSave(i,j)
            end select

            xy_SurfRH(i,j) =                                    &
              & - xyr_RadSFlux(i,j,0)                           &
              & - xyr_RadLFlux(i,j,0)                           &
              & - xyr_HeatFlux(i,j,0)                           &
              & - xy_SurfLatentHeatFlux(i,j)                    &
!              & + xy_DeepSubSurfHeatFlux(i,j)                   &
              & + xyr_SoilHeatFlux(i,j,0)                       &
              & - xy_LatHeatFluxByOtherSpc(i,j)                 &
              & - LatentHeatLocal * DelSurfSol / ( 2.0_DP * DelTime )

          end if
        end do
      end do
      !
      xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
      xya_SoilTempVec = xya_ArgSoilTempVec


      ! 温度の計算
      ! Calculate temperature and specific humidity
      !
      do l = -1, 1
        do k = 1, kslmax
          xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
        end do
        k = 0
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
        do k = 1, kmax
          xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
        end do
      end do
      !
      call PhyImplLUDecomp3(              &
        & xyaa_TempSoilTempLUMtx,         & ! (inout)
        & imax * jmax, kmax + 1 + kslmax  & ! (in)
        )
      !
      do k = 1, kslmax
        xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
      end do
      k = 0
      xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
      do k = 1, kmax
        xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
      end do
      !
      call PhyImplLUSolve3(                   &
        & xya_DelTempSoilTempLUVec,           & ! (inout)
        & xyaa_TempSoilTempLUMtx,             & ! (in)
        & 1, imax * jmax , kmax + 1 + kslmax  & ! (in)
        & )

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

            if ( xy_FlagCalc(i,j) ) then

              select case ( xy_IndexCalcMethod(i,j) )
              case ( IndexLand )
                xyz_DSoilTempDt(i,j,k) = &
                  & xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
              case default
                xyz_DSoilTempDt(i,j,k) = 0.0_DP
              end select

            end if

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

          if ( xy_FlagCalc(i,j) ) then

            select case ( xy_IndexCalcMethod(i,j) )
            case ( IndexLand )
              ! land
              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSeaIce )
!!$              ! sea ice
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSlabOcean )
!!$              ! slab ocean
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexOceanPresSST )
!!$              ! open ocean
!!$              xy_DSurfTempDt(i,j) = 0.0_DP
            case default
              call MessageNotify( 'E', module_name, 'Unexpected Error.' )
            end select

          end if

        end do
      end do
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
            end if
          end do
        end do
      end do
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime

            select case ( xy_IndexMeltOrFreeze(i,j) )
            case ( IndexMelt )
              ! all ice/snow melt
              DelSurfSol =   xy_SurfSolATentativeSave(i,j)
            case ( IndexFreeze )
              ! all soil moisture freeze (= negative melt of ice/snow)
              DelSurfSol = - xy_SurfLiqATentativeSave(i,j)
            end select

            xy_LatHeatFluxBySnowMelt(i,j) = &
              & LatentHeatLocal * DelSurfSol / ( 2.0_DP * DelTime )
            xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
              & + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
          end if
        end do
      end do


!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xy_FlagCalc(i,j) ) then
!!$            if ( xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime ) < xy_TempCond(i,j) ) then
!!$              call MessageNotify( 'M', module_name, &
!!$                & 'Surface temperature is lower than condensation temperature, %f < %f.', &
!!$                & d = (/ xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime ), xy_TempCond(i,j) /) )
!!$            end if
!!$          end if
!!$        end do
!!$      end do


      !----------


      ! Calculation for a land point with prescribed temperature
      !
      !----------
      ! Surface fluxes used below
      !----------
!!$      xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
!!$      xy_SurfRadLFlux = xyr_RadLFlux(:,:,0)                                           &
!!$        &   + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime )        &
!!$        &   + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
!!$      xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0)                                &
!!$        & - xyr_SoilTempTransCoef(:,:,0)                                               &
!!$        &   * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0d0 * DelTime )
!!$      xy_SurfSensHeatFlux =                                                       &
!!$        & xyr_HeatFlux(:,:,0)                                                     &
!!$        &   - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef                     &
!!$        &     * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1)                           &
!!$        &       - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
!!$      !
!!$      xy_LatHeatFluxBySnowMelt =           &
!!$        & - xy_SurfRadSFlux                &
!!$        & - xy_SurfRadLFlux                &
!!$        & - xy_SurfSensHeatFlux            &
!!$        & - xy_SurfLatentHeatFlux          &
!!$        & + xy_SurfSoilHeatCondFlux        &
!!$        & - xy_LatHeatFluxByOtherSpc       &
!!$        & - xy_SurfHeatCapacity * xy_DSurfTempDt
      !
      do j = 1, jmax
        do i = 0, imax-1
          select case ( xy_IndexCalcMethod(i,j) )
          case ( IndexLandWithPresTs )
            SurfSolATentative = xy_SurfSolATentativeSave(i,j)
            LatHeatFluxBySnowMelt =           &
              & - xy_SurfRadSFlux(i,j)                &
              & - xy_SurfRadLFlux(i,j)                &
              & - xy_SurfSensHeatFlux(i,j)            &
              & - xy_SurfLatentHeatFlux(i,j)          &
!!$              & + xy_SurfSoilHeatCondFlux(i,j)        & ! Comment out since the flux is not calculated correctly.
              & - xy_LatHeatFluxByOtherSpc(i,j)       &
              & - xy_SurfHeatCapacity(i,j) * xy_DSurfTempDt(i,j)

            if ( SurfSolATentative > 0.0_DP ) then
              ! Ice exists on the ground.
              !   Calculation is performed only when freezing and melting
              xy_FlagCalc(i,j) = .true.
              xy_LatHeatFluxBySnowMelt(i,j) = LatHeatFluxBySnowMelt
!!$            else if ( xy_LatHeatFluxBySnowMelt(i,j) < 0.0_DP ) then
            else if ( LatHeatFluxBySnowMelt < 0.0_DP ) then
              ! Ice does not exist on the ground.
              !   Calculation is performed only when freezing
              xy_FlagCalc(i,j) = .true.
              xy_LatHeatFluxBySnowMelt(i,j) = LatHeatFluxBySnowMelt
            else
              xy_FlagCalc(i,j) = .false.
            end if
          case default
            xy_FlagCalc(i,j) = .false.
          end select
        end do
      end do
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
            SurfSolATentative = xy_SurfSolATentativeSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal * ( 2.0_DP * DelTime )
            if ( SurfSolATentative < 0.0_DP ) then
              xy_LatHeatFluxBySnowMelt(i,j) = &
                & LatentHeatLocal * SurfSolATentative / ( 2.0_DP * DelTime )
            end if
!!$              xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
!!$                & - SurfSolATentative / ( 2.0_DP * DelTime )
!!$              xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
!!$                & + SurfSolATentative / ( 2.0_DP * DelTime )
            xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
              & + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
          end if
        end do
      end do


    end if




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

  end subroutine OLD_PhyImplSDHV3IceSnowPhase

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

  subroutine PhyImplSDHV3IceSnowPhaseChgCor(         &
    & IndexSpc,                                      & ! (in)
    & xy_Ps,                                         & ! (in)
    & xyr_HeatFlux, xy_SurfLatentHeatFlux,           & ! (in)
    & xyr_SoilHeatFlux,                              & ! (in)
    & xyr_SoilTempTransCoef,                         & ! (in)
    & xyr_RadSFlux, xyr_RadLFlux,                    & ! (in)
    & xy_DeepSubSurfHeatFlux,                        & ! (in)
    & xy_SurfTemp, xyz_SoilTemp,                     & ! (in)
    & xy_SurfLiqB, xy_SurfSolB,                      & ! (in)
    & xy_SurfHeatCapacity,                           & ! (in)
    & xy_SoilHeatCap, xy_SoilHeatDiffCoef,           & ! (in)
    & xy_IndexCalcMethod,                            & ! (in)
    & xyra_DelRadLFlux,                              & ! (in)
    & xyz_Exner, xyr_Exner,                          & ! (in)
    & xy_SurfTempTransCoef,                          & ! (in)
    & xy_LatHeatFluxByOtherSpc,                      & ! (in)
    & xyza_ArgTempMtx, xyz_ArgTempVec,               & ! (in)
    & xyaa_ArgSurfMtx, xy_ArgSurfRH,                 & ! (in)
    & xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec,       & ! (in)
    & xyz_DTempDt,                                   & ! (in)
    & xy_DSurfTempDt,                                & ! (inout)
    & xyz_DSoilTempDt,                               & ! (inout)
    & xy_DSurfLiqDt, xy_DSurfSolDt,                  & ! (inout)
    & xy_LatHeatFluxBySnowMelt                       & ! (out)
    & )
    !
    ! 融雪による時間変化率の修正を行います. 
    !
    ! Correction of tendencies due to melt of snow. 
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_SSDepth,   &    ! subsurface grid on interface of layer
      & z_SSDepth         ! subsurface grid at midpoint of layer

    ! 物理定数設定
    ! 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
      & LatentHeatFusion
                              ! $ L $ [J kg-1] . 
                              ! 融解の潜熱. 
                              ! Latent heat of fusion

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

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


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

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

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


    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: IndexSpc
    real(DP), intent(in):: xy_Ps(0:imax-1, 1:jmax)
                              ! 
                              ! Surface pressure

    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux

    real(DP), intent(in):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface

    real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
                              ! 土壌中の熱フラックス (W m-2)
                              ! Heat flux in sub-surface soil (W m-2)
    real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
                              ! 輸送係数：土壌温度.
                              ! Transfer coefficient: soil temperature


    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux

    real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
                              ! 地中熱フラックス. 
                              ! "Deep subsurface heat flux"
                              ! Heat flux at the bottom of surface/soil layer.
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xyz_SoilTemp (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度 (K)
                              ! Soil temperature (K)

    real(DP), intent(in):: xy_SurfLiqB (0:imax-1, 1:jmax)
                              ! 
                              ! Surface liquid amount
    real(DP), intent(in):: xy_SurfSolB (0:imax-1, 1:jmax)
                              ! 積雪量.
                              ! Surface snow amount.

    real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
                              ! 地表熱容量. 
                              ! Surface heat capacity
    real(DP), intent(in ):: xy_SoilHeatCap (0:imax-1, 1:jmax)
                              ! 土壌熱容量 (J K-1 kg-1)
                              ! Specific heat of soil (J K-1 kg-1)
    real(DP), intent(in ):: xy_SoilHeatDiffCoef (0:imax-1, 1:jmax)
                              ! 土壌熱伝導係数 (J m-3 K-1)
                              ! Heat conduction coefficient of soil (J m-3 K-1)

    integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method

    real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature

    real(DP), intent(in   ):: xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP), intent(in   ):: xyz_ArgTempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP), intent(in   ):: xyaa_ArgSurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance
    real(DP), intent(in   ):: xy_ArgSurfRH(0:imax-1,1:jmax)


    real(DP), intent(in   ):: xyaa_ArgSoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP), intent(in   ):: xya_ArgSoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature

    real(DP), intent(in   ):: xy_LatHeatFluxByOtherSpc(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by other specie

    real(DP), intent(inout):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(inout):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)
    real(DP), intent(inout):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
                              ! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
                              ! Temperature tendency (K s-1)

    real(DP), intent(inout):: xy_DSurfLiqDt (0:imax-1, 1:jmax)
                              ! 土壌温度時間変化率 (kg m-2 s-1)
                              ! Soil temperature tendency (kg m-2 s-1)
    real(DP), intent(inout):: xy_DSurfSolDt (0:imax-1, 1:jmax)
                              ! 積雪率時間変化率 (kg m-2 s-1)
                              ! Surface snow amount tendency (kg m-2 s-1)
    real(DP), intent(out  ):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by melt (variable only for debug)

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

    real(DP):: xy_TempCond          (0:imax-1, 1:jmax)
    real(DP):: xy_MajCompPressSatB  (0:imax-1, 1:jmax)
    real(DP):: xy_DMajCompPressSatDT(0:imax-1, 1:jmax)

    logical :: xy_FlagCalc(0:imax-1, 1:jmax)
    integer :: xy_IndexMeltOrFreeze(0:imax-1, 1:jmax)
    integer, parameter :: IndexOthers = 0
    integer, parameter :: IndexMelt   = 1
    integer, parameter :: IndexFreeze = 2

    real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance
    real(DP):: xy_SurfRH(0:imax-1,1:jmax)


    real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature

    real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
                              ! $ T, Tg $ の時間変化.
                              ! Tendency of $ T $ and $ Tg |

    real(DP):: LatentHeatLocal
    real(DP):: LatentHeatFluxByMelt
    real(DP):: SenHeatFluxA
    real(DP):: LatHeatFluxA
    real(DP):: CondHeatFluxA
    real(DP):: ValueAlpha
    real(DP):: SurfTempATentative
    real(DP):: SoilTempATentative
    real(DP):: SurfLiqATentative
    real(DP):: xy_SurfLiqATentativeSave(0:imax-1, 1:jmax)
    real(DP):: SurfSolATentative
    real(DP):: xy_SurfSolATentativeSave(0:imax-1, 1:jmax)
    real(DP):: DelSurfSol

    real(DP) :: xy_TempMajCompCond(0:imax-1, 1:jmax)
    real(DP) :: SurfMajCompIceATentative

    real(DP) :: xy_SurfRadSFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfRadLFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSensHeatFlux    (0:imax-1, 1:jmax)
    real(DP) :: xy_SeaIceHeatCondFlux  (0:imax-1, 1:jmax)
    real(DP) :: xy_HeatingTendency     (0:imax-1, 1:jmax)

    real(DP) :: LatHeatFluxBySnowMelt

    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


    ! 実行文 ; Executable statement
    !

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

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


    !
    ! check flag of snow melt
    !
    if ( ( .not. FlagBucketModel ) .or. ( .not. FlagSnow ) ) then
      xy_LatHeatFluxBySnowMelt = 0.0_DP
      return
    end if


    if ( kslmax == 0 ) then

      call MessageNotify( 'E', module_name, &
        & 'kslmax <= 0 in PhyImplSDHV3IceSnowPhaseChgCor.' )

    else

      xy_DSurfLiqDtSave = xy_DSurfLiqDt
      xy_DSurfSolDtSave = xy_DSurfSolDt


      select case ( IndexSpc )
      case ( IndexSpcMajComp )
        call SaturateMajorCompCondTemp( &
          & xy_Ps,                      & ! (in)
          & xy_TempCond                 & ! (out)
          & )
        call SaturateMajorCompPressSat(     &
          & xy_SurfTemp,                    &  ! (in)
          & xy_MajCompPressSatB             &  ! (out)
          & )
        call SaturateMajorCompDPressSatDT(  &
          & xy_SurfTemp,                    &  ! (in)
          & xy_DMajCompPressSatDT           &  ! (out)
          & )
        LatentHeatLocal = SaturateMajorCompInqLatentHeat()
      case ( IndexSpcH2O )
        xy_TempCond     = TempCondWater
        LatentHeatLocal = LatentHeatFusion
      case default
        call MessageNotify( 'E', module_name, &
          & 'Undefined IndexSpc, %d.', i = (/ IndexSpc /) )
      end select


      xy_SurfLiqATentativeSave = xy_SurfLiqB + xy_DSurfLiqDt * ( 2.0_DP * DelTime )
      xy_SurfSolATentativeSave = xy_SurfSolB + xy_DSurfSolDt * ( 2.0_DP * DelTime )

      !----------
      ! A case that a part of snow/ice melt or soil moisture freeze
      !----------

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

          select case ( xy_IndexCalcMethod(i,j) )
          case ( IndexLand )

            SurfTempATentative = xy_SurfTemp(i,j)          &
              & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
            SurfSolATentative = xy_SurfSolATentativeSave(i,j)

            if ( &
              & ( SurfSolATentative  > 0.0_DP           ) .and. &
              & ( SurfTempATentative > xy_TempCond(i,j) )       &
              & ) then
              xy_FlagCalc         (i,j) = .true.
            else
              xy_FlagCalc         (i,j) = .false.
            end if

          case default
            xy_FlagCalc         (i,j) = .false.
          end select

        end do
      end do

      ! Freeze
      select case ( IndexSpc )
      case ( IndexSpcMajComp )

        do j = 1, jmax
          do i = 0, imax-1
            select case ( xy_IndexCalcMethod(i,j) )
            case ( IndexLand )
              SurfTempATentative = xy_SurfTemp(i,j)          &
                & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
              if ( SurfTempATentative < xy_TempCond(i,j) ) then
                xy_FlagCalc         (i,j) = .true.
              end if
            end select
          end do
        end do

!!$      case ( IndexSpcH2O )
!!$
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            select case ( xy_IndexCalcMethod(i,j) )
!!$            case ( IndexLand )
!!$              SurfTempATentative = xy_SurfTemp(i,j)          &
!!$                & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
!!$              SurfLiqATentative = xy_SurfLiqATentativeSave(i,j)
!!$              if ( &
!!$                & ( SurfLiqATentative  > 0.0_DP           ) .and. &
!!$                & ( SurfTempATentative < xy_TempCond(i,j) )       &
!!$                & ) then
!!$                xy_FlagCalc         (i,j) = .true.
!!$              end if
!!$            end select
!!$
!!$          end do
!!$        end do
!!$
      end select


      xyza_TempMtx     = xyza_ArgTempMtx
      xyz_TempVec = xyz_ArgTempVec
      !
      xyaa_SurfMtx     = xyaa_ArgSurfMtx
      xy_SurfRH        = xy_ArgSurfRH
      !
      select case ( IndexSpc )
      case ( IndexSpcMajComp )
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              xyaa_SurfMtx(i,j,0, 0) = xyaa_SurfMtx(i,j,0, 0) &
                & + LatentHeatLocal / Grav * xy_DMajCompPressSatDT(i,j) &
                &   / ( 2.0_DP * DelTime )
              xy_SurfRH   (i,j)      = xy_SurfRH(i,j)                                  &
                & - LatentHeatLocal / Grav * ( xy_MajCompPressSatB(i,j) - xy_Ps(i,j) ) &
                &   / ( 2.0_DP * DelTime )
            end if
          end do
        end do
      case ( IndexSpcH2O )
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
              xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
              xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
              xy_SurfRH   (i,j)      = xy_TempCond(i,j) - xy_SurfTemp(i,j)
            end if
          end do
        end do
      end select
      !
      xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
      xya_SoilTempVec = xya_ArgSoilTempVec


      ! 温度の計算
      ! Calculate temperature and specific humidity
      !
      do l = -1, 1
        do k = 1, kslmax
          xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
        end do
        k = 0
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
        do k = 1, kmax
          xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
        end do
      end do

      call PhyImplLUDecomp3(              &
        & xyaa_TempSoilTempLUMtx,         & ! (inout)
        & imax * jmax, kmax + 1 + kslmax  & ! (in)
        )

      do k = 1, kslmax
        xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
      end do
      k = 0
      xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
      do k = 1, kmax
        xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
      end do

      call PhyImplLUSolve3(                   &
        & xya_DelTempSoilTempLUVec,           & ! (inout)
        & xyaa_TempSoilTempLUMtx,             & ! (in)
        & 1, imax * jmax , kmax + 1 + kslmax  & ! (in)
        & )


      do k = 1, kslmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              select case ( xy_IndexCalcMethod(i,j) )
              case ( IndexLand )
                xyz_DSoilTempDt(i,j,k) = &
                  & xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
              case default
                xyz_DSoilTempDt(i,j,k) = 0.0_DP
              end select
            end if
          end do
        end do
      end do
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
            select case ( xy_IndexCalcMethod(i,j) )
            case ( IndexLand )
              ! land
              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSeaIce )
!!$              ! sea ice
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSlabOcean )
!!$              ! slab ocean
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexOceanPresSST )
!!$              ! open ocean
!!$              xy_DSurfTempDt(i,j) = 0.0_DP
            case default
              call MessageNotify( 'E', module_name, 'Unexpected Error.' )
            end select
          end if
        end do
      end do
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
            end if
          end do
        end do
      end do


      !----------
      ! Surface fluxes used below
      !----------
      xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
      xy_SurfRadLFlux = xyr_RadLFlux(:,:,0)                                           &
        &   + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime )        &
        &   + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
      xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0)                                &
        & - xyr_SoilTempTransCoef(:,:,0)                                               &
        &   * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0d0 * DelTime )
      xy_SurfSensHeatFlux =                                                       &
        & xyr_HeatFlux(:,:,0)                                                     &
        &   - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef                     &
        &     * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1)                           &
        &       - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )

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

          if ( xy_FlagCalc(i,j) ) then

            xy_LatHeatFluxBySnowMelt(i,j) =           &
              & - xy_SurfRadSFlux(i,j)                &
              & - xy_SurfRadLFlux(i,j)                &
              & - xy_SurfSensHeatFlux(i,j)            &
              & - xy_SurfLatentHeatFlux(i,j)          &
              & + xy_SurfSoilHeatCondFlux(i,j)        &
              & - xy_LatHeatFluxByOtherSpc(i,j)       &
              & - xy_SurfHeatCapacity(i,j) * xy_DSurfTempDt(i,j)

            xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
              & + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal

!!$            if ( xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * ( 2.0_DP * DelTime ) < 0.0_DP ) then
!!$              call MessageNotify( 'M', module_name, &
!!$                & 'Surface snow amount is negative %f, %f.', &
!!$                & d = (/ xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * ( 2.0_DP * DelTime ), xy_SurfSnowB(i,j) /) )
!!$            end if

          else
            xy_LatHeatFluxBySnowMelt(i,j) = 0.0_DP
          end if

        end do
      end do

      !----------
      ! A case that all snow melt or soil moisture freeze
      !----------

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

          select case ( xy_IndexCalcMethod(i,j) )
          case ( IndexLand )

            if ( xy_FlagCalc(i,j) ) then
              SurfSolATentative = xy_SurfSolB(i,j) &
                & + xy_DSurfSolDt(i,j) * 2.0_DP * DelTime
              if ( SurfSolATentative < 0.0_DP ) then
                xy_FlagCalc(i,j) = .true.
                xy_IndexMeltOrFreeze(i,j) = IndexMelt
              else
                xy_FlagCalc(i,j) = .false.
                xy_IndexMeltOrFreeze(i,j) = IndexOthers
              end if

!!$              select case ( IndexSpc )
!!$              case ( IndexSpcH2O )
!!$                SurfLiqATentative = xy_SurfLiqB(i,j) &
!!$                  & + xy_DSurfLiqDt(i,j) * 2.0_DP * DelTime
!!$                if ( SurfLiqATentative < 0.0_DP ) then
!!$                  xy_FlagCalc(i,j) = .true.
!!$                  xy_IndexMeltOrFreeze(i,j) = IndexFreeze
!!$                end if
!!$              end select

            else
              xy_FlagCalc(i,j) = .false.
              xy_IndexMeltOrFreeze(i,j) = IndexOthers
            end if

          case default
            xy_FlagCalc(i,j) = .false.
            xy_IndexMeltOrFreeze(i,j) = IndexOthers
          end select

        end do
      end do


      xyza_TempMtx     = xyza_ArgTempMtx
      xyz_TempVec      = xyz_ArgTempVec
      !
      xyaa_SurfMtx     = xyaa_ArgSurfMtx
      xy_SurfRH        = xy_ArgSurfRH
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
            select case ( xy_IndexMeltOrFreeze(i,j) )
            case ( IndexMelt )
              ! all ice/snow melt
              DelSurfSol =   xy_SurfSolATentativeSave(i,j)
            case ( IndexFreeze )
              ! all soil moisture freeze (= negative melt of ice/snow)
              DelSurfSol = - xy_SurfLiqATentativeSave(i,j)
            end select

            xy_SurfRH(i,j) =                                    &
              & - xyr_RadSFlux(i,j,0)                           &
              & - xyr_RadLFlux(i,j,0)                           &
              & - xyr_HeatFlux(i,j,0)                           &
              & - xy_SurfLatentHeatFlux(i,j)                    &
!              & + xy_DeepSubSurfHeatFlux(i,j)                   &
              & + xyr_SoilHeatFlux(i,j,0)                       &
              & - xy_LatHeatFluxByOtherSpc(i,j)                 &
              & - LatentHeatLocal * DelSurfSol / ( 2.0_DP * DelTime )

          end if
        end do
      end do
      !
      xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
      xya_SoilTempVec = xya_ArgSoilTempVec


      ! 温度の計算
      ! Calculate temperature and specific humidity
      !
      do l = -1, 1
        do k = 1, kslmax
          xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
        end do
        k = 0
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
        do k = 1, kmax
          xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
        end do
      end do
      !
      call PhyImplLUDecomp3(              &
        & xyaa_TempSoilTempLUMtx,         & ! (inout)
        & imax * jmax, kmax + 1 + kslmax  & ! (in)
        )
      !
      do k = 1, kslmax
        xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
      end do
      k = 0
      xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
      do k = 1, kmax
        xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
      end do
      !
      call PhyImplLUSolve3(                   &
        & xya_DelTempSoilTempLUVec,           & ! (inout)
        & xyaa_TempSoilTempLUMtx,             & ! (in)
        & 1, imax * jmax , kmax + 1 + kslmax  & ! (in)
        & )

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

            if ( xy_FlagCalc(i,j) ) then

              select case ( xy_IndexCalcMethod(i,j) )
              case ( IndexLand )
                xyz_DSoilTempDt(i,j,k) = &
                  & xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
              case default
                xyz_DSoilTempDt(i,j,k) = 0.0_DP
              end select

            end if

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

          if ( xy_FlagCalc(i,j) ) then

            select case ( xy_IndexCalcMethod(i,j) )
            case ( IndexLand )
              ! land
              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSeaIce )
!!$              ! sea ice
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexSlabOcean )
!!$              ! slab ocean
!!$              xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$            case ( IndexOceanPresSST )
!!$              ! open ocean
!!$              xy_DSurfTempDt(i,j) = 0.0_DP
            case default
              call MessageNotify( 'E', module_name, 'Unexpected Error.' )
            end select

          end if

        end do
      end do
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_FlagCalc(i,j) ) then
              xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
            end if
          end do
        end do
      end do
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime

            select case ( xy_IndexMeltOrFreeze(i,j) )
            case ( IndexMelt )
              ! all ice/snow melt
              DelSurfSol =   xy_SurfSolATentativeSave(i,j)
            case ( IndexFreeze )
              ! all soil moisture freeze (= negative melt of ice/snow)
              DelSurfSol = - xy_SurfLiqATentativeSave(i,j)
            end select

            xy_LatHeatFluxBySnowMelt(i,j) = &
              & LatentHeatLocal * DelSurfSol / ( 2.0_DP * DelTime )
            xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
              & + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
          end if
        end do
      end do


!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xy_FlagCalc(i,j) ) then
!!$            if ( xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime ) < xy_TempCond(i,j) ) then
!!$              call MessageNotify( 'M', module_name, &
!!$                & 'Surface temperature is lower than condensation temperature, %f < %f.', &
!!$                & d = (/ xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime ), xy_TempCond(i,j) /) )
!!$            end if
!!$          end if
!!$        end do
!!$      end do


      !----------


      ! Calculation for a land point with prescribed temperature
      !
      !----------
      ! Surface fluxes used below
      !----------
!!$      xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
!!$      xy_SurfRadLFlux = xyr_RadLFlux(:,:,0)                                           &
!!$        &   + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime )        &
!!$        &   + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
!!$      xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0)                                &
!!$        & - xyr_SoilTempTransCoef(:,:,0)                                               &
!!$        &   * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0d0 * DelTime )
!!$      xy_SurfSensHeatFlux =                                                       &
!!$        & xyr_HeatFlux(:,:,0)                                                     &
!!$        &   - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef                     &
!!$        &     * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1)                           &
!!$        &       - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
!!$      !
!!$      xy_LatHeatFluxBySnowMelt =           &
!!$        & - xy_SurfRadSFlux                &
!!$        & - xy_SurfRadLFlux                &
!!$        & - xy_SurfSensHeatFlux            &
!!$        & - xy_SurfLatentHeatFlux          &
!!$        & + xy_SurfSoilHeatCondFlux        &
!!$        & - xy_LatHeatFluxByOtherSpc       &
!!$        & - xy_SurfHeatCapacity * xy_DSurfTempDt
      !
      do j = 1, jmax
        do i = 0, imax-1
          select case ( xy_IndexCalcMethod(i,j) )
          case ( IndexLandWithPresTs )
            SurfSolATentative = xy_SurfSolATentativeSave(i,j)
            LatHeatFluxBySnowMelt =           &
              & - xy_SurfRadSFlux(i,j)                &
              & - xy_SurfRadLFlux(i,j)                &
              & - xy_SurfSensHeatFlux(i,j)            &
              & - xy_SurfLatentHeatFlux(i,j)          &
!!$              & + xy_SurfSoilHeatCondFlux(i,j)        & ! Comment out since the flux is not calculated correctly.
              & - xy_LatHeatFluxByOtherSpc(i,j)       &
              & - xy_SurfHeatCapacity(i,j) * xy_DSurfTempDt(i,j)

            if ( SurfSolATentative > 0.0_DP ) then
              ! Ice exists on the ground.
              !   Calculation is performed only when freezing and melting
              xy_FlagCalc(i,j) = .true.
              xy_LatHeatFluxBySnowMelt(i,j) = LatHeatFluxBySnowMelt
!!$            else if ( xy_LatHeatFluxBySnowMelt(i,j) < 0.0_DP ) then
            else if ( LatHeatFluxBySnowMelt < 0.0_DP ) then
              ! Ice does not exist on the ground.
              !   Calculation is performed only when freezing
              xy_FlagCalc(i,j) = .true.
              xy_LatHeatFluxBySnowMelt(i,j) = LatHeatFluxBySnowMelt
            else
              xy_FlagCalc(i,j) = .false.
            end if
          case default
            xy_FlagCalc(i,j) = .false.
          end select
        end do
      end do
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
            SurfSolATentative = xy_SurfSolATentativeSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal * ( 2.0_DP * DelTime )
            if ( SurfSolATentative < 0.0_DP ) then
              xy_LatHeatFluxBySnowMelt(i,j) = &
                & LatentHeatLocal * SurfSolATentative / ( 2.0_DP * DelTime )
            end if
!!$              xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
!!$                & - SurfSolATentative / ( 2.0_DP * DelTime )
!!$              xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
!!$                & + SurfSolATentative / ( 2.0_DP * DelTime )
            xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
              & + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
          end if
        end do
      end do


    end if




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

  end subroutine PhyImplSDHV3IceSnowPhaseChgCor

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

  subroutine PhyImplSDHV3SeaIceCorrection(                                 &
    & xy_IndexCalcMethod,                                                  & ! (in)
    & xy_SeaIceThickness,                                                  & ! (in)
    & xyz_Exner, xyr_Exner,                                                & ! (in)
    & xy_SurfTemp,                                                         & ! (in)
    & xyr_HeatFlux,                                                        & ! (in)
    & xyr_SoilHeatFlux,                                                    & ! (in)
    & xy_SurfTempTransCoef, xyr_SoilTempTransCoef,                         & ! (in)
    & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux,                        & ! (in)
    & xy_SurfLatentHeatFlux,                                               & ! (in)
    & xyza_ArgTempMtx, xyz_ArgTempVec,                                     & ! (in)
    & xyaa_ArgSurfMtx, xy_ArgSurfRH,                                       & ! (in)
    & xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec,                             & ! (in)
    & xyz_DTempDt,                                                         & ! (inout)
    & xy_DSurfTempDt,                                                      & ! (inout)
    & xyz_DSoilTempDt,                                                     & ! (inout)
    & xy_LatHeatFluxBySeaIceMelt                                           & ! (out)
    & )
    !
    ! 融雪による時間変化率の修正を行います. 
    !
    ! Correction of tendencies due to melt of snow. 
    !

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

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & CpDry
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure

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

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in   ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
                              ! 
                              ! Index for calculation method
    real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              !
                              ! Sea ice thickness
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in   ) :: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature

    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
                              ! 土壌中の熱フラックス (W m-2)
                              ! Heat flux in sub-surface soil (W m-2)

    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
                              ! 輸送係数：土壌温度.
                              ! Transfer coefficient: soil temperature

    real(DP), intent(in   ):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in   ):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(in   ):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave

    real(DP), intent(in   ):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface

    real(DP), intent(in   ):: xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP), intent(in   ):: xyz_ArgTempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP), intent(in   ):: xyaa_ArgSurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance
    real(DP), intent(in   ):: xy_ArgSurfRH(0:imax-1,1:jmax)


    real(DP), intent(in   ):: xyaa_ArgSoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP), intent(in   ):: xya_ArgSoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature

    real(DP), intent(inout):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(inout) :: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度変化率 (K s-1)
                              ! Surface temperature tendency (K s-1)
    real(DP), intent(inout):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
                              ! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
                              ! Temperature tendency (K s-1)
    real(DP), intent(out  ) :: xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax)


    ! 作業変数
    ! Work variables
    !

    real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 温度陰解ベクトル. 
                              ! Implicit vector about temperature
    real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
                              ! 惑星表面エネルギー収支用陰解行列
                              ! Implicit matrix for surface energy balance
    real(DP):: xy_SurfRH(0:imax-1,1:jmax)


    real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
                              ! 土壌温度拡散方程式の行列
                              ! Matrix for diffusion equation of soil temperature
    real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
                              ! 土壌温度拡散方程式のベクトル
                              ! Vector for diffusion equation of soil temperature
    real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
                              ! $ T, Tg $ の時間変化.
                              ! Tendency of $ T $ and $ Tg |


    logical :: xy_FlagSeaIceMelt(0:imax-1, 1:jmax)

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

    real(DP) :: SurfTempATentative

    real(DP) :: xy_SurfRadSFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfRadLFlux        (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
    real(DP) :: xy_SurfSensHeatFlux    (0:imax-1, 1:jmax)
    real(DP) :: xy_SeaIceHeatCondFlux  (0:imax-1, 1:jmax)
    real(DP) :: xy_HeatingTendency     (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


    ! 実行文 ; Executable statement
    !

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


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


    !
    ! check flag of snow melt
    !
    if ( ( .not. FlagBucketModel ) .or. ( .not. FlagSnow ) ) then
      xy_LatHeatFluxBySeaIceMelt = 0.0_DP
      return
    end if


    xy_TempCond = TempCondWater

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

        SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
        if ( &
          & ( xy_IndexCalcMethod(i,j) == IndexSeaIce ) .and.   &
          & ( SurfTempATentative > xy_TempCond(i,j) )                  &
          & ) then
          xy_FlagSeaIceMelt(i,j) = .true.
        else
          xy_FlagSeaIceMelt(i,j) = .false.
        end if

      end do
    end do


    xyza_TempMtx     = xyza_ArgTempMtx
    xyz_TempVec      = xyz_ArgTempVec
    !
    xyaa_SurfMtx     = xyaa_ArgSurfMtx
    xy_SurfRH        = xy_ArgSurfRH
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_FlagSeaIceMelt(i,j) ) then
          xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
          xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
          xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
          xy_SurfRH   (i,j)      = xy_TempCond(i,j) - xy_SurfTemp(i,j)
        end if
      end do
    end do
    !
    xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
    xya_SoilTempVec = xya_ArgSoilTempVec



    ! 温度の計算
    ! Calculate temperature and specific humidity
    !
    do l = -1, 1
      do k = 1, kslmax
        xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
      end do
      k = 0
      xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
      do k = 1, kmax
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
      end do
    end do

    call PhyImplLUDecomp3(              &
      & xyaa_TempSoilTempLUMtx,         & ! (inout)
      & imax * jmax, kmax + 1 + kslmax  & ! (in)
      )

    do k = 1, kslmax
      xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
    end do
    k = 0
    xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
    do k = 1, kmax
      xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
    end do

    call PhyImplLUSolve3(                   &
      & xya_DelTempSoilTempLUVec,           & ! (inout)
      & xyaa_TempSoilTempLUMtx,             & ! (in)
      & 1, imax * jmax , kmax + 1 + kslmax  & ! (in)
      & )



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

          if ( xy_FlagSeaIceMelt(i,j) ) then

            select case ( xy_IndexCalcMethod(i,j) )
!!$            case ( IndexLand )
!!$              xyz_DSoilTempDt(i,j,k) = &
!!$                & xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
            case default
              xyz_DSoilTempDt(i,j,k) = 0.0_DP
            end select

          end if

        end do
      end do
    end do
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_FlagSeaIceMelt(i,j) ) then
          select case ( xy_IndexCalcMethod(i,j) )
!!$          case ( IndexLand )
!!$            ! land
!!$            xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
          case ( IndexSeaIce )
            ! sea ice
            xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$          case ( IndexSlabOcean )
!!$            ! slab ocean
!!$            xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$          case ( IndexOceanPresSST )
!!$            ! open ocean
!!$            xy_DSurfTempDt(i,j) = 0.0_DP
          case default
            call MessageNotify( 'E', module_name, 'Unexpected Error.' )
          end select
        end if
      end do
    end do
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1

          if ( xy_FlagSeaIceMelt(i,j) ) then
            xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
          end if

        end do
      end do
    end do


    xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
    xy_SurfRadLFlux = xyr_RadLFlux(:,:,0)                                           &
      &   + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime )        &
      &   + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
    xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0)                                &
      & - xyr_SoilTempTransCoef(:,:,0)                                               &
      &   * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0_DP * DelTime )
    xy_SurfSensHeatFlux =                                                       &
      & xyr_HeatFlux(:,:,0)                                                     &
      &   - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef                     &
      &     * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1)                           &
      &       - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
    xy_SeaIceHeatCondFlux =                                                           &
      & - SeaIceThermCondCoef                                                         &
      &   * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) &
      &   / xy_SeaIceThickness
    xy_HeatingTendency = SeaIceVolHeatCap * xy_SeaIceThickness * xy_DSurfTempDt

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

        if ( xy_FlagSeaIceMelt(i,j) ) then
          xy_LatHeatFluxBySeaIceMelt(i,j) =          &
            & - (                                    &
            &     (                                  &
            &         xy_SurfRadSFlux(i,j)           &
            &       + xy_SurfRadLFlux(i,j)           &
            &       + xy_SurfSensHeatFlux(i,j)       &
            &       + xy_SurfLatentHeatFlux(i,j)     &
            &     )                                  &
            &     - xy_SeaIceHeatCondFlux(i,j)       &
            &   )                                    &
            & - xy_HeatingTendency(i,j)
        else
          xy_LatHeatFluxBySeaIceMelt(i,j) = 0.0_DP
        end if

      end do
    end do


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


  end subroutine PhyImplSDHV3SeaIceCorrection

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

  subroutine PhyImplSDHV3Init(                                      &
    & ArgFlagBucketModel, ArgFlagSnow,                              & ! (in)
    & ArgFlagSlabOcean, ArgFlagMajCompPhaseChange, CondMajCompName  & ! (in)
    & )

    !
    ! phy_implicit モジュールの初期化を行います. 
    ! NAMELIST#phy_implicit_nml の読み込みはこの手続きで行われます. 
    !
    ! "phy_implicit" module is initialized. 
    ! "NAMELIST#phy_implicit_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

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

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

    ! 地表面フラックスユーティリティ
    ! Surface flux utility routines
    !
    use surface_flux_util, only : SurfaceFluxUtilInit


    ! 宣言文 ; Declaration statements
    !

    logical     , intent(in ) :: ArgFlagBucketModel
                              ! flag for use of bucket model
    logical     , intent(in ) :: ArgFlagSnow
                              ! flag for treating snow
    logical     , intent(in ) :: ArgFlagSlabOcean
                              ! flag for use of slab ocean
    logical     , intent(in ) :: ArgFlagMajCompPhaseChange
                              ! flag for use of major component phase change
    character(*), intent(in ) :: CondMajCompName


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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /phy_implicit_sdh_V3_nml/ &
      & SOHeatCapacity, &          ! Slab ocean heat capacity (J m-2 K-1)
      & NumMaxItr,      &          ! Number of interation
      & TempItrCrit,    &
      & FlagSublimation
          !
          ! デフォルト値については初期化手続 "phy_implicit#PhyImplInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "phy_implicit#PhyImplInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( phy_implicit_sdh_V3_inited ) return


    ! Set flag for bucket model
    FlagBucketModel = ArgFlagBucketModel

    ! Set flag for treating snow
    FlagSnow = ArgFlagSnow

    ! Set flag for slab ocean
    FlagSlabOcean = ArgFlagSlabOcean

    ! Set flag for major component phase change
    FlagMajCompPhaseChange = ArgFlagMajCompPhaseChange

    ! デフォルト値の設定
    ! Default values settings
    !

    SOHeatCapacity = 4.187e3_DP * 1.0e3_DP * 60.0_DP
                         ! 4.187d3 (J (kg K)-1) * 1.0d3 (kg m-3) * 60.0d0 (m)

    NumMaxItr   = 50
    TempItrCrit = 0.05_DP

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

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


    ! Initialization of modules used in this model
    !

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

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

    ! 地表面フラックスユーティリティ
    ! Surface flux utility routines
    !
    call SurfaceFluxUtilInit


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  SOHeatCapacity  = %f', d = (/ SOHeatCapacity /) )
    call MessageNotify( 'M', module_name, '  NumMaxItr       = %d', i = (/ NumMaxItr /) )
    call MessageNotify( 'M', module_name, '  TempItrCrit     = %f', d = (/ TempItrCrit /) )
    call MessageNotify( 'M', module_name, '  FlagSublimation = %b', l = (/ FlagSublimation /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    phy_implicit_sdh_V3_inited = .true.

  end subroutine PhyImplSDHV3Init

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

end module phy_implicit_sdh_V3