!= 地下熱伝導モデルを用いた場合の陰解法による時間積分
!
!= 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.f90,v 1.15 2014/05/07 09:39:19 murashin Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2009. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module phy_implicit_sdh
  !
  != 地下熱伝導モデルを用いた場合の陰解法による時間積分
  !
  != 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_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 :: PhyImplSDHTendency
  public :: PhyImplSDHSetMethodFromMatthews
  public :: PhyImplSDHCorSOTempBySnowMelt
  public :: PhyImplSDHInit


  ! 公開変数
  ! 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

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



  ! Index for calculation method
  !
  integer, parameter :: IndexOceanPresSST = 10   ! Ocean with prescribed SST
  integer, parameter :: IndexSlabOcean    = 11   ! Slab ocean
  integer, parameter :: IndexSeaIce       = 12   ! Sea ice
  integer, parameter :: IndexLand         = 13   ! Land

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

  character(*), parameter:: module_name = 'phy_implicit_sdh'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: phy_implicit_sdh.f90,v 1.15 2014/05/07 09:39:19 murashin Exp $'
                              ! モジュールのバージョン
                              ! Module version


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


contains

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

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

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

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


    ! 宣言文 ; Declaration statements
    !

    integer , intent(in ) :: xy_SurfCond       (0:imax-1, 1:jmax)
                              ! 地表状態. 
                              ! Surface condition
    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_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) >= 1 ) then
          ! land
          xy_IndexCalcMethod(i,j) = IndexLand
        else
          if ( xy_SeaIceConc(i,j) > SeaIceThreshold ) then
            ! sea ice
            xy_IndexCalcMethod(i,j) = IndexSeaIce
          else if ( FlagSlabOcean ) then
            ! slab ocean
            xy_IndexCalcMethod(i,j) = IndexSlabOcean
          else
            ! open ocean
            xy_IndexCalcMethod(i,j) = IndexOceanPresSST
          end if
        end if
      end do
    end do


  end subroutine PhyImplSDHSetMethodFromMatthews

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

  subroutine PhyImplSDHCorSOTempBySnowMelt(  &
    & 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
    !

    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_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
        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 * SeaIceThickness )
        case ( IndexSlabOcean )
          xy_SurfTemp(i,j) = xy_SurfTemp(i,j)                              &
            & + LatentHeatFusion * xy_SurfSnowFlux(i,j) * 2.0_DP * DelTime &
            &   / SOHeatCapacity
        case ( IndexOceanPresSST )
        case default
          call MessageNotify( 'E', module_name, 'This index is inappropriate.' )
        end select
      end do
    end do


  end subroutine PhyImplSDHCorSOTempBySnowMelt

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

  subroutine PhyImplSDHTendency(                               &
    & xy_IndexCalcMethod,                                      & ! (in)
    & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in)
    & 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_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
    !

    ! 座標データ設定
    ! 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, &
      & SeaIceThreshold,    &
      & SeaIceThickness,    &
      & 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):: 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(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_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):: 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):: xy_SurfQVapSat (0:imax-1, 1:jmax)
!!$                              ! 地表飽和比湿. 
!!$                              ! Saturated specific humidity on surface
!!$    real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
!!$                              ! 地表飽和比湿変化. 
!!$                              ! Saturated specific humidity tendency on surface

    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)

    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_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 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
      do k = 1, kmax
        xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
      end do
    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 * SeaIceThickness / ( 2.0_DP * DelTime )  &
            & + CpDry * xy_SurfTempTransCoef(i,j)                         &
            & + xyra_DelRadLFlux(i,j,0,0)                                 &
            & + SeaIceThermCondCoef / SeaIceThickness
          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 ( IndexOceanPresSST )
          ! open ocean
          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 ) / SeaIceThickness
        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 ( IndexOceanPresSST )
          ! open ocean
          xy_SurfRH(i,j) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    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.0_DP * DelTime )
      xyz_DVDt(:,:,k) = xyz_DVDt(:,:,k) / ( 2.0_DP * DelTime )
    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 ( IndexOceanPresSST )
          ! open ocean
          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

              xy_DSurfSnowDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
              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_DSoilMoistDt(i,j) = 0.0_DP
              else
                xy_DSurfSnowDt (i,j) = 0.0_DP
                xy_DSoilMoistDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
              end if

            end if

          end do
        end do
      else
        ! Evaporation is subtracted from soil moisture
        !
        xy_DSoilMoistDt = - xyrf_QMixFlux(:,:,0,IndexH2OVap)
        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

!!$    call PhyImplSDHSnowMeltCorrection(                 &
!!$      & 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_SurfMajCompIceB,                            & ! (in)
!!$      & xy_SurfHeatCapacity,                           & ! (in)
!!$      & xy_SoilHeatCap, xy_SoilHeatDiffCoef,           & ! (in)
!!$      & xy_IndexCalcMethod,                            & ! (in)
!!$      & xyra_DelRadLFlux,                              & ! (in)
!!$      & xyz_Exner, xyr_Exner,                          & ! (in)
!!$      & xy_SurfTempTransCoef,                          & ! (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_DSurfMajCompIceDt,                          & ! (inout)
!!$      & xy_DSoilMoistDt, xy_DSurfSnowDt,               & ! (inout)
!!$      & xy_LatHeatFluxByMajCompIceSubl,                & ! (out)
!!$      & xy_LatHeatFluxBySnowMelt                       & ! (out)
!!$      & )


    if ( FlagMajCompPhaseChange ) then

      xy_DAtmMassDt        = 0.0_DP
      xy_DSurfMajCompIceDt = 0.0_DP
      ! A dummy value
      !
      xy_LatHeatFluxByOtherSpc = 0.0_DP

      call PhyImplSDHIceSnowPhaseChangeCor(              &
        & 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_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

    call PhyImplSDHIceSnowPhaseChangeCor(              &
      & 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_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)
      & )



    call PhyImplSDHSeaIceCorrection(                               &
      & xy_IndexCalcMethod,                                        & ! (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 PhyImplSDHChkConservation(                            &
      & xy_IndexCalcMethod,                                    &
      & xyr_Press, xyz_Exner, xyr_Exner,                       &
      & xy_SurfTemp,                                           &
      & xy_SurfHeatCapacity, xy_SoilHeatCap,                   &
      & xyr_HeatFlux, xyrf_QMixFlux, 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 PhyImplSDHTendency

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

  subroutine PhyImplSDHChkConservation(                       &
    & xy_IndexCalcMethod,                                     &
    & xyr_Press, xyz_Exner, xyr_Exner,                        &
    & xy_SurfTemp,                                            &
    & xy_SurfHeatCapacity, xy_SoilHeatCap,                    &
    & xyr_HeatFlux, xyrf_QMixFlux, 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, &
      & SeaIceThreshold,    &
      & SeaIceThickness,    &
      & TempBelowSeaIce

    ! 宣言文 ; Declaration statements
    !

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

    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_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_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 ) &
      &   / 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_Residual =                        &
      &   xy_SurfRadSFlux                &
      & + xy_SurfRadLFlux                &
      & + xy_SurfSensHeatFlux            &
      & + xy_SurfLatentHeatFlux          &
      & - xy_SurfSoilHeatCondFlux        &
      & + xy_LatHeatFluxByMajCompIceSubl &
      & + xy_LatHeatFluxBySnowMelt
    !
    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, &
        & '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.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 ( 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, &
        & 'Slab ocean heating res.    : %f.', d = (/ MaxResidual /) )
    end if


    !-----
    ! Sea ice heating
    !
    xy_SumAtmRate = SeaIceVolHeatCap * 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 ( IndexOceanPresSST )
          ! open ocean
        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) ) &
      & - 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.0d-10 ) 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 ( IndexOceanPresSST )
          ! open ocean
        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 PhyImplSDHChkConservation

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

  subroutine PhyImplSDHIceSnowPhaseChangeCor(        &
    & 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_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_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)

    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):: SurfSolATentative
    real(DP):: xy_SurfSolATentativeSave(0:imax-1, 1:jmax)

    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)

    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_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 PhyImplSDHSnowMeltCorrection.' )

    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_SurfSolATentativeSave = xy_SurfSolB + xy_DSurfSolDt * ( 2.0_DP * DelTime )


      !----------
      ! A case that a part of snow/ice melts
      !----------

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

          if ( xy_IndexCalcMethod(i,j) == IndexLand ) then

            SurfTempATentative = xy_SurfTemp(i,j)          &
              & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * 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

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

        end do
      end do


      select case ( IndexSpc )
      case ( IndexSpcMajComp )

        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
              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 if
          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



      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 )

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

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

          if ( xy_IndexCalcMethod(i,j) == IndexLand ) then

            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.
              else
                xy_FlagCalc(i,j) = .false.
              end if
            else
              xy_FlagCalc(i,j) = .false.
            end if

          else
            xy_FlagCalc(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_FlagCalc(i,j) ) then
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
            SurfSolATentative = xy_SurfSolATentativeSave(i,j)
            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 * SurfSolATentative / ( 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
            SurfSolATentative = xy_SurfSolATentativeSave(i,j)

            xy_LatHeatFluxBySnowMelt(i,j) = &
              & LatentHeatLocal * 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


      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


      !----------

    end if



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

  end subroutine PhyImplSDHIceSnowPhaseChangeCor

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

  subroutine PhyImplSDHSnowMeltCorrection(           &
    & 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_SurfSnowB,                                  & ! (in)
    & 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)
    & 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_DSurfMajCompIceDt,                          & ! (inout)
    & xy_DSoilMoistDt, xy_DSurfSnowDt,               & ! (inout)
    & xy_LatHeatFluxByMajCompIceSubl,                & ! (out)
    & 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


    ! 宣言文 ; Declaration statements
    !

    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_SurfSnowB (0:imax-1, 1:jmax)
                              ! 積雪量.
                              ! Surface snow amount.

    real(DP), intent(in):: xy_SurfMajCompIceB(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice 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(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_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Major component ice amount tendency (kg m-2 s-1)

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

    real(DP), intent(out  ):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
                              !
                              ! Latent heat flux by major component ice sublimation
                              ! (variable only for debug)
    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_DSoilMoistDtSave(0:imax-1, 1:jmax)
    real(DP):: xy_DSurfSnowDtSave (0:imax-1, 1:jmax)

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

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

    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):: SurfSnowATentative
    real(DP):: xy_SurfSnowATentativeSave(0:imax-1, 1:jmax)

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


    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_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_LatHeatFluxByMajCompIceSubl = 0.0_DP
      xy_LatHeatFluxBySnowMelt       = 0.0_DP
      return
    end if


    if ( kslmax == 0 ) then

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

    else

      xy_DSoilMoistDtSave = xy_DSoilMoistDt
      xy_DSurfSnowDtSave  = xy_DSurfSnowDt


      xy_LatHeatFluxByMajCompIceSubl = 0.0_DP


      xy_TempCond     = TempCondWater
      LatentHeatLocal = LatentHeatFusion


      xy_SurfSnowATentativeSave = xy_SurfSnowB + xy_DSurfSnowDt * ( 2.0_DP * DelTime )


      !----------
      ! A case that a part of snow melts
      !----------

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

          if ( xy_IndexCalcMethod(i,j) == IndexLand ) then

            SurfTempATentative = xy_SurfTemp(i,j)          &
              & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
            SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)

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

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

        end do
      end do


      xyza_TempMtx     = xyza_ArgTempMtx
      xyz_TempVec = xyz_ArgTempVec
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
            k = 1
            xyza_TempMtx(i,j,k,-1) = 0.0_DP
            xyz_TempVec(i,j,k) = - ( xyr_HeatFlux(i,j,k) - xyr_HeatFlux(i,j,k-1) )
            xyz_TempVec(i,j,k) = xyz_TempVec(i,j,k) &
              & + CpDry * xy_SurfTempTransCoef(i,j) &
              &   * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
          end if
        end do
      end do
      !
      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, 0) = 1.0_DP
            xy_SurfRH(i,j) =                                      &
              & - xyr_RadSFlux(i,j,0)                             &
              & - xyr_RadLFlux(i,j,0)                             &
              & - xyr_HeatFlux(i,j,0)                             &
              & - xy_SurfLatentHeatFlux(i,j)                      &
              & + xyr_SoilHeatFlux(i,j,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)                )  &
              &   * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
          end if
        end do
      end do
      !
      xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
      xya_SoilTempVec = xya_ArgSoilTempVec
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagCalc(i,j) ) then
            k = 1
            xyaa_SoilTempMtx(i,j,k,-1) = 0.0_DP
            xya_SoilTempVec(i,j,k) = - ( xyr_SoilHeatFlux(i,j,k) - xyr_SoilHeatFlux(i,j,k-1) )
            xya_SoilTempVec(i,j,k) = xya_SoilTempVec(i,j,k) &
              & + xyr_SoilTempTransCoef(i,j,k-1) &
              &    * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
          end if
        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

            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
            xy_DSurfTempDt(i,j) = &
              & ( xy_TempCond(i,j) - xy_SurfTemp(i,j) ) / ( 2.0_DP * DelTime )
          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
            xy_LatHeatFluxBySnowMelt(i,j) = xya_DelTempSoilTempLUVec(i,j,0)

!!$            xy_DSurfSnowDt(i,j) = - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
!!$            xy_DSoilMoistDt(i,j) = - xy_DSurfSnowDt(i,j)

            xy_DSurfSnowDt(i,j) = xy_DSurfSnowDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSoilMoistDt(i,j) = xy_DSoilMoistDtSave(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 melts
      !----------

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

          if ( xy_IndexCalcMethod(i,j) == IndexLand ) then

            if ( xy_FlagCalc(i,j) ) then
              SurfSnowATentative = xy_SurfSnowB(i,j) &
                & + xy_DSurfSnowDt(i,j) * 2.0_DP * DelTime
!!$              SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
!!$              if ( SurfSnowATentative < 0.0_DP ) then
!!$              SurfTempATentative = &
!!$                & xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime )
!!$              if ( SurfTempATentative < xy_TempCond(i,j) ) then
              if ( SurfSnowATentative < 0.0_DP ) then
                xy_FlagCalc(i,j) = .true.
              else
                xy_FlagCalc(i,j) = .false.
              end if
            else
              xy_FlagCalc(i,j) = .false.
            end if

          else
            xy_FlagCalc(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_FlagCalc(i,j) ) then
!!$            SurfSnowATentative = xy_SurfSnowB(i,j)         &
!!$              & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
            SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
            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)                       &
              & - LatentHeatLocal * SurfSnowATentative / ( 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.
            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
            SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)

            xy_LatHeatFluxBySnowMelt(i,j) = &
              & LatentHeatLocal * SurfSnowATentative / ( 2.0_DP * DelTime )
            xy_DSurfSnowDt(i,j) = xy_DSurfSnowDtSave(i,j) &
              & - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
            xy_DSoilMoistDt(i,j) = xy_DSoilMoistDtSave(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


      !----------

    end if



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

  end subroutine PhyImplSDHSnowMeltCorrection

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

  subroutine PhyImplSDHSeaIceCorrection(                                   &
    & xy_IndexCalcMethod,                                                  & ! (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, &
      & SeaIceThreshold    , &
      & SeaIceVolHeatCap   , &
      & SeaIceThickness    , &
      & 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):: 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_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. * 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 ) &
      &   / SeaIceThickness
    xy_HeatingTendency = SeaIceVolHeatCap * 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 PhyImplSDHSeaIceCorrection

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

  subroutine PhyImplSDHInit(                                        &
    & 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

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


    ! 宣言文 ; 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_nml/ &
      SOHeatCapacity          ! Slab ocean heat capacity (J m-2 K-1)
          !
          ! デフォルト値については初期化手続 "phy_implicit#PhyImplInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "phy_implicit#PhyImplInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( phy_implicit_sdh_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)


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

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


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


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


    phy_implicit_sdh_inited = .true.

  end subroutine PhyImplSDHInit

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

end module phy_implicit_sdh
