subroutine PhyImplSDHV3CorSOTempBySnowMelt( 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
    !
    real(DP) :: xy_SeaIceThickness(0:imax-1, 1:jmax)

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


    ! 実行文 ; Executable statement
    !

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


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

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


    !
    ! Set sea ice thickness
    !
    xy_SeaIceThickness = SeaIceThickness

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


  end subroutine PhyImplSDHV3CorSOTempBySnowMelt

  subroutine PhyImplSDHV3Init( 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

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

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

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


    ! 宣言文 ; Declaration statements
    !

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


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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /phy_implicit_sdh_V3_nml/ SOHeatCapacity, NumMaxItr, TempItrCrit, FlagSublimation
          !
          ! デフォルト値については初期化手続 "phy_implicit#PhyImplInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "phy_implicit#PhyImplInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( phy_implicit_sdh_V3_inited ) return


    ! Set flag for bucket model
    FlagBucketModel = ArgFlagBucketModel

    ! Set flag for treating snow
    FlagSnow = ArgFlagSnow

    ! Set flag for slab ocean
    FlagSlabOcean = ArgFlagSlabOcean

    ! Set flag for major component phase change
    FlagMajCompPhaseChange = ArgFlagMajCompPhaseChange

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

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

    NumMaxItr   = 50
    TempItrCrit = 0.05_DP

    FlagSublimation = .false.

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = phy_implicit_sdh_V3_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

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


    ! Initialization of modules used in this model
    !

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

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

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


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


    phy_implicit_sdh_V3_inited = .true.

  end subroutine PhyImplSDHV3Init

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

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

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


    ! 宣言文 ; Declaration statements
    !

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


    ! 作業変数
    ! Work variables
    !

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


    ! 実行文 ; Executable statement
    !

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


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

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


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


  end subroutine PhyImplSDHV3SetMethodMatthews

  subroutine PhyImplSDHV3Tendency( xy_IndexCalcMethod, xy_BucketFlagOceanGrid, xy_SnowFrac, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xyz_TempB, xy_SurfTemp, xyz_SoilTemp, xyzf_QMixB, xy_SurfHumidCoef, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyr_SoilTempTransCoef, xy_SurfMajCompIceB, xy_SoilMoistB, 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
    !

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

    ! 座標データ設定
    ! 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, LatentHeat

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

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

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

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

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

    ! 宣言文 ; Declaration statements
    !

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

    integer:: iitr

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: l               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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


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

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

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

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

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

    !
    ! Set sea ice thickness
    !
    xy_SeaIceThickness = SeaIceThickness


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


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

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

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

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


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


    ! Calculation for thermal diffusion
    !

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


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

    xyz_TempSave    = xyz_TempB
    xy_SurfTempSave = xy_SurfTemp

    ! iteration
    iitr = 1
    loop_itr : do

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


      ! Tendencies of atmospheric, surface, and subsurface temperatures, and 
      ! atmospheric water vapor are solved with a fixed surface water vapor flux. 
      !
      n = IndexH2OVap
      xy_SurfH2OVapFlux = xyrf_QMixFlux(:,:,0,n) - xy_SurfHumidCoef * xy_SurfQVapTransCoef * ( xyzf_DQMixDt(:,:,1,n) - xy_SurfDQVapSatDTemp * xy_DSurfTempDt ) * 2.0_DP * DelTime
      ! Limit surface flux not to be negative atmospheric content
      ! IMPORTANT : Now, only the water vapor flux is restricted.
      call SurfaceFluxUtilLimitFlux( ( 2.0_DP * DelTime ), xyzf_QMixB, xyr_Press, xy_SurfH2OVapFlux )
      ! Calculation of latent heat flux
      xy_SurfLatentHeatFlux = LatentHeat * xy_SurfH2OVapFlux
      !
      if ( FlagBucketModel ) then
        ! バケツモデルのための地表面フラックス修正
        ! Modification of surface flux for bucket model
        !
        call BucketModEvapAndLatentHeatFlux( xy_BucketFlagOceanGrid, xy_SoilMoistB, xy_SurfSnowB, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux )
      end if
      !
      call PhyImplSDHV3TendencyHeatCore( xy_IndexCalcMethod, xy_SeaIceThickness, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFlux, 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_SoilMoistB, xy_SurfSnowB, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DPsDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt )


      ! Check

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

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

      a_LocalMax(1) = MaxTempInc
      a_LocalMax(2) = MaxSurfTempInc
      call MPIWrapperFindMaxVal( 2, a_LocalMax, a_GlobalMax )
      MaxTempInc     = a_GlobalMax(1)
      MaxSurfTempInc = a_GlobalMax(2)

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

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

      xyz_TempSave    = xyz_TempA
      xy_SurfTempSave = xy_SurfTempA

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

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


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


  end subroutine PhyImplSDHV3Tendency

  subroutine OLD_PhyImplSDHV3IceSnowPhase( IndexSpc, xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfLiqB, 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 : SaturateMajorCompCondTemp, SaturateMajorCompInqLatentHeat


    ! 宣言文 ; Declaration statements
    !

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

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

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

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


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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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


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

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

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

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

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

    real(DP) :: LatHeatFluxBySnowMelt

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


    ! 実行文 ; Executable statement
    !

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

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


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


    if ( kslmax == 0 ) then

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

    else

      xy_DSurfLiqDtSave = xy_DSurfLiqDt
      xy_DSurfSolDtSave = xy_DSurfSolDt


      select case ( IndexSpc )
      case ( IndexSpcMajComp )
        call SaturateMajorCompCondTemp( xy_Ps, 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_SurfLiqATentativeSave = xy_SurfLiqB + xy_DSurfLiqDt * ( 2.0_DP * DelTime )
      xy_SurfSolATentativeSave = xy_SurfSolB + xy_DSurfSolDt * ( 2.0_DP * DelTime )

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

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

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

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

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

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

        end do
      end do

      ! Freeze
      select case ( IndexSpc )
      case ( IndexSpcMajComp )

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

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


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


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

      call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, 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


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

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

          if ( xy_FlagCalc(i,j) ) then

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

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

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

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

        end do
      end do

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

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

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

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

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

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

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

        end do
      end do


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

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

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


      ! 温度の計算
      ! Calculate temperature and specific humidity
      !
      do l = -1, 1
        do k = 1, kslmax
          xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
        end do
        k = 0
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
        do k = 1, kmax
          xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
        end do
      end do
      !
      call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, 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

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

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


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


      !----------


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

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


    end if




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

  end subroutine OLD_PhyImplSDHV3IceSnowPhase

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

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: r_SSDepth, 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, TempBelowSeaIce

    ! 宣言文 ; Declaration statements
    !

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

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

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

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

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

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

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

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

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

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


    ! 作業変数
    ! Work variables
    !

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

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

    real(DP) :: MaxResidual

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


    ! 実行文 ; Executable statement
    !

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


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


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

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

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


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


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


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


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


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


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


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


  end subroutine PhyImplSDHV3ChkConservation

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

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: r_SSDepth, 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, LatentHeat

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

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

    ! 宣言文 ; Declaration statements
    !

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

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

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

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

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

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


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

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

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

    real(DP) :: MaxResidual

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


    ! 実行文 ; Executable statement
    !

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


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


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

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


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

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

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


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


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


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


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


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


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


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


  end subroutine PhyImplSDHV3ChkConservationTQ

  subroutine PhyImplSDHV3IceSnowPhaseChgCor( IndexSpc, xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfLiqB, 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: Grav, 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 : SaturateMajorCompCondTemp, SaturateMajorCompPressSat, SaturateMajorCompDPressSatDT, SaturateMajorCompInqLatentHeat


    ! 宣言文 ; Declaration statements
    !

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

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

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

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


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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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


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

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

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

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

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

    real(DP) :: LatHeatFluxBySnowMelt

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


    ! 実行文 ; Executable statement
    !

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

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


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


    if ( kslmax == 0 ) then

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

    else

      xy_DSurfLiqDtSave = xy_DSurfLiqDt
      xy_DSurfSolDtSave = xy_DSurfSolDt


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


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

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

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

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

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

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

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

        end do
      end do

      ! Freeze
      select case ( IndexSpc )
      case ( IndexSpcMajComp )

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

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


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


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

      call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, 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


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

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

          if ( xy_FlagCalc(i,j) ) then

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

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

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

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

        end do
      end do

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

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

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

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

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

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

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

        end do
      end do


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

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

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


      ! 温度の計算
      ! Calculate temperature and specific humidity
      !
      do l = -1, 1
        do k = 1, kslmax
          xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
        end do
        k = 0
        xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
        do k = 1, kmax
          xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
        end do
      end do
      !
      call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, 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

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

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


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


      !----------


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

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


    end if




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

  end subroutine PhyImplSDHV3IceSnowPhaseChgCor

  subroutine PhyImplSDHV3SeaIceCorrection( xy_IndexCalcMethod, xy_SeaIceThickness, 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, SeaIceVolHeatCap   , TempBelowSeaIce

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

    ! 宣言文 ; Declaration statements
    !

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

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

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

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

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

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

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


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

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


    ! 作業変数
    ! Work variables
    !

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


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


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

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

    real(DP) :: SurfTempATentative

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


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


    ! 実行文 ; Executable statement
    !

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


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


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


    xy_TempCond = TempCondWater

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

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

      end do
    end do


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



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

    call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, 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.0_DP * DelTime )
          end if

        end do
      end do
    end do


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

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

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

      end do
    end do


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


  end subroutine PhyImplSDHV3SeaIceCorrection

  subroutine PhyImplSDHV3TendencyHeatCore( xy_IndexCalcMethod, xy_SeaIceThickness, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFlux, 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_SoilMoistB, xy_SurfSnowB, 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, 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):: xy_SeaIceThickness(0:imax-1, 1:jmax)
                              !
                              ! Sea ice thickness

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: l               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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



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

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



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

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

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

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

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


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

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

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

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

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



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

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

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

    end if

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


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

      end do
    end do

    do j = 1, jmax
      do i = 0, imax-1
        select case ( xy_IndexCalcMethod(i,j) )
        case ( IndexLand )
          ! land
          xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - xy_SurfLatentHeatFlux(i,j) + 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 ) / xy_SeaIceThickness(i,j)
        case ( IndexSlabOcean )
          ! slab ocean
          xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - xy_SurfLatentHeatFlux(i,j)                       !&
!              & + xy_DeepSubSurfHeatFlux(i,j)
        case ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xy_SurfRH(i,j) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do


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

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

    call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, 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 ( IndexPresTs, IndexLandWithPresTs )
          ! prescribed surface temperature
          xy_DSurfTempDt(i,j) = 0.0_DP
        case default
          call MessageNotify( 'E', module_name, 'Unexpected Error.' )
        end select
      end do
    end do
    do k = 1, kmax
      xyz_DTempDt(:,:,k) = xya_DelTempSoilTempLUVec(:,:,k) / ( 2.0_DP * DelTime )
    end do


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

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

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

            else

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

            end if

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


    ! Temporarily set
    !
    xy_DSurfMajCompIceDt = 0.0_DP


    if ( FlagMajCompPhaseChange ) then

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

      call PhyImplSDHV3IceSnowPhaseChgCor( IndexSpcMajComp, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfMajCompLiqB, 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

    if ( FlagSublimation ) then
      ! If sublimation is considered, the melt of snow/ice is not calculated.
      xy_LatHeatFluxBySnowMelt = 0.0_DP
    else
      ! Else, the melt of snow/ice is calculated.
      call PhyImplSDHV3IceSnowPhaseChgCor( IndexSpcH2O, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SoilMoistB, 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 )
    end if



    call PhyImplSDHV3SeaIceCorrection( xy_IndexCalcMethod, xy_SeaIceThickness, 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 PhyImplSDHV3ChkConservation( xy_IndexCalcMethod, xy_SeaIceThickness, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHeatCapacity, xy_SoilHeatCap, xyr_HeatFlux, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_LatHeatFluxByMajCompIceSubl, xy_LatHeatFluxBySnowMelt, xy_LatHeatFluxBySeaIceMelt, xy_DeepSubSurfHeatFlux, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_DPsDt, xy_DSurfMajCompIceDt )


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

  end subroutine PhyImplSDHV3TendencyHeatCore

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, GasRDry, LatentHeat

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

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

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

    ! 宣言文 ; Declaration statements
    !

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

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

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

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

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

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

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

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

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

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


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

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

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: l               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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



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

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


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


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

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

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

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

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


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

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

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

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

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



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

      end do
    end do

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


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

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

    call PhyImplLUDecomp3( xyaa_TempQVapLUMtx, imax * jmax, kmax + 1 + kmax )

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

    call PhyImplLUSolve3( xya_DelTempQVapLUVec, xyaa_TempQVapLUMtx, 1, imax * jmax , kmax + 1 + kmax )

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


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


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

  end subroutine PhyImplSDHV3TendencyHeatTQCore

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, GasRDry
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air

    ! 時刻管理
    ! 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
    !

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

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

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

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


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

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


    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: l               ! 行列用 DO ループ用作業変数
                              ! Work variables for DO loop of matrices
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

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


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



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

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



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

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

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

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

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

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

    call PhyImplLUDecomp3( xyza_UVLUMtx, 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. * DelTime )
      xyz_DVDt(:,:,k) = xyz_DVDt(:,:,k) / ( 2. * DelTime )
    end do


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

  end subroutine PhyImplSDHV3TendencyMomCore