!= 蝨ｰ荳狗�莨晏ｰ弱Δ繝�Ν繧堤畑縺�◆蝣ｴ蜷医�髯ｰ隗｣豕輔↓繧医ｋ譎る俣遨榊�
!
!= 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