subroutine PhyImplSDHCorSOTempBySnowMelt( xy_IndexCalcMethod, xy_SurfSnowFlux, xy_SurfTemp )
    !
    !
    !
    ! 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 PhyImplSDHInit( ArgFlagBucketModel, ArgFlagSnow, ArgFlagSlabOcean, ArgFlagMajCompPhaseChange, CondMajCompName )

    !
    ! 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, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = phy_implicit_sdh_nml, 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 )
    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

  subroutine PhyImplSDHSetMethodFromMatthews( xy_SurfCond, xy_SeaIceConc, xy_IndexCalcMethod )
    !
    !
    !
    ! 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 PhyImplSDHTendency( xy_IndexCalcMethod, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfHumidCoef, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyr_SoilTempTransCoef, xy_SurfMajCompIceB, xy_SurfSnowB, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DPsDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt )
    !
    ! 時間変化率の計算を行います. 
    !
    ! Calculate tendencies. 
    !

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, 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, TimeN, 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) + 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) - 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, 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, xyza_UVLUMtx, 1, imax * jmax, kmax ) ! (in)

    call PhyImplLUSolve3( xyz_DVDt, xyza_UVLUMtx, 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, imax * jmax, kmax + 1 + kslmax )

    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, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )

    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, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfMajCompIceB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xy_LatHeatFluxByOtherSpc, xyza_TempMtx, xyz_TempVec, xyaa_SurfMtx, xy_SurfRH, xyaa_SoilTempMtx, xya_SoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DAtmMassDt, xy_DSurfMajCompIceDt, xy_LatHeatFluxByMajCompIceSubl )

    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, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfSnowB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xy_LatHeatFluxByOtherSpc, xyza_TempMtx, xyz_TempVec, xyaa_SurfMtx, xy_SurfRH, xyaa_SoilTempMtx, xya_SoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_LatHeatFluxBySnowMelt )



    call PhyImplSDHSeaIceCorrection( xy_IndexCalcMethod, xyz_Exner, xyr_Exner, xy_SurfTemp, xyr_HeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_SurfLatentHeatFlux, xyza_TempMtx, xyz_TempVec, xyaa_SurfMtx, xy_SurfRH, xyaa_SoilTempMtx, xya_SoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_LatHeatFluxBySeaIceMelt )


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

    call PhyImplLUDecomp3( xyza_QMixLUMtx, imax * jmax, kmax )

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

      call PhyImplLUSolve3( xyz_DelQMixLUVec, xyza_QMixLUMtx, 1, imax * jmax , kmax )

      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, z_SSDepth         ! subsurface grid at midpoint of layer

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, 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, xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfSolB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xy_LatHeatFluxByOtherSpc, xyza_ArgTempMtx, xyz_ArgTempVec, xyaa_ArgSurfMtx, xy_ArgSurfRH, xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSurfLiqDt, xy_DSurfSolDt, xy_LatHeatFluxBySnowMelt )
    !
    ! 融雪による時間変化率の修正を行います. 
    !
    ! Correction of tendencies due to melt of snow. 
    !

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: CpDry, 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, TimeN, 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 : SaturateMajorCompCalcCondTemp, 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 SaturateMajorCompCalcCondTemp( xy_Ps, xy_TempCond )
        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, imax * jmax, kmax + 1 + kslmax )

      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, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )


      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) + 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, imax * jmax, kmax + 1 + kslmax )
      !
      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, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )

      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 PhyImplSDHSeaIceCorrection( xy_IndexCalcMethod, xyz_Exner, xyr_Exner, xy_SurfTemp, xyr_HeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_SurfLatentHeatFlux, xyza_ArgTempMtx, xyz_ArgTempVec, xyaa_ArgSurfMtx, xy_ArgSurfRH, xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_LatHeatFluxBySeaIceMelt )
    !
    ! 融雪による時間変化率の修正を行います. 
    !
    ! Correction of tendencies due to melt of snow. 
    !

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

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, 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, imax * jmax, kmax + 1 + kslmax )

    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, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )



    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 PhyImplSDHSnowMeltCorrection( xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfSnowB, xy_SurfMajCompIceB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xyza_ArgTempMtx, xyz_ArgTempVec, xyaa_ArgSurfMtx, xy_ArgSurfRH, xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_LatHeatFluxByMajCompIceSubl, xy_LatHeatFluxBySnowMelt )
    !
    ! 融雪による時間変化率の修正を行います. 
    !
    ! Correction of tendencies due to melt of snow. 
    !

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: CpDry, 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, TimeN, 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, imax * jmax, kmax + 1 + kslmax )

      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, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )


      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) + 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, imax * jmax, kmax + 1 + kslmax )
      !
      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, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )

      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