subroutine PhyImplicitClose( phy_impl, err )
    !
    ! PHYIMPL 型の変数の終了処理を行います.
    ! なお, 与えられた *phy_impl* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Deconstructor of "PHYIMPL".
    ! Note that if *phy_impl* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYIMPL), intent(inout):: phy_impl
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyImplicitClose'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_impl % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYIMPL'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  "PHYIMPL" の設定の消去
    !  Clear the settings for "PHYIMPL"
    !-----------------------------------------------------------------
    deallocate( phy_impl % xy_SurfHeatCapacity )
    deallocate( phy_impl % xy_SurfCondition )
    deallocate( phy_impl % xy_GroundTempFlux )

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
    phy_impl % initialized = .false.
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyImplicitClose

  subroutine PhyImplicitCreate( phy_impl, imax, jmax, kmax, Grav, Cp, EL, DelTime, xy_SurfHeatCapacity, xy_SurfCondition, xy_GroundTempFlux, nmlfile, err )
    !
    ! PHYIMPL 型の変数の初期設定を行います.
    ! 他のサブルーチンを使用する前に必ずこのサブルーチンによって
    ! PHYIMPL 型の変数を初期設定してください.
    !
    ! なお, 与えられた *phy_impl* が既に初期設定されている場合,
    ! プログラムはエラーを発生させます.
    !
    ! NAMELIST を利用する場合には引数 *nmlfile* に NAMELIST ファイル名
    ! を与えてください. NAMELIST 変数群の詳細に関しては 
    ! NAMELIST#phy_implicit_nml を参照してください. 
    !
    ! Constructor of "PHYIMPL".
    ! Initialize *phy_impl* by this subroutine, 
    ! before other procedures are used, 
    !
    ! Note that if *phy_impl* is already initialized 
    ! by this procedure, error is occurred.
    !
    ! In order to use NAMELIST, specify a NAMELIST filename to 
    ! argument *nmlfile*. See "NAMELIST#phy_implicit_nml"
    ! for details about a NAMELIST group.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_present, only: present_and_not_empty, present_and_true
    use dc_message, only: MessageNotify
    use dc_error, only: StoreError, DC_NOERR, DC_EALREADYINIT, DC_EARGLACK, DC_ENEGATIVE, DC_ENOFILEREAD
    implicit none
    type(PHYIMPL), intent(inout):: phy_impl
    integer, intent(in):: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer, intent(in):: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude
    integer, intent(in):: kmax ! 鉛直層数. 
                   ! Number of vertical level
    real(DP), intent(in):: Grav      ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP), intent(in):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
    real(DP), intent(in):: EL        ! $ L $ .      水の凝結の潜熱. Latent heat of condensation of water vapor
    real(DP), intent(in):: DelTime    ! $ \Delta t $ . タイムステップ. Time step
    real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 0:jmax-1)
                              ! 地表熱容量. 
                              ! Surface heat capacity
    integer, intent(in):: xy_SurfCondition (0:imax-1, 0:jmax-1)
                              ! 地表状態. 
                              ! Surface condition
    real(DP), intent(in):: xy_GroundTempFlux (0:imax-1, 0:jmax-1)
                              ! 地中熱フラックス. 
                              ! Ground temperature flux
    character(*), intent(in), optional:: nmlfile
                              ! NAMELIST ファイルの名称. 
                              ! この引数に空文字以外を与えた場合, 
                              ! 指定されたファイルから 
                              ! NAMELIST 変数群を読み込みます. 
                              ! ファイルを読み込めない場合にはエラーを
                              ! 生じます.
                              !
                              ! NAMELIST 変数群の詳細に関しては 
                              ! NAMELIST#phy_implicit_nml 
                              ! を参照してください. 
                              !
                              ! NAMELIST file name. 
                              ! If nonnull character is specified to
                              ! this argument, 
                              ! NAMELIST group name is loaded from the 
                              ! file. 
                              ! If the file can not be read, 
                              ! an error occurs.
                              ! 
                              ! See "NAMELIST#phy_implicit_nml" 
                              ! for details about a NAMELIST group.
                              ! 
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyImplicitCreate'
  continue
    call BeginSub( subname, version )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( phy_impl % initialized ) then
      stat = DC_EALREADYINIT
      cause_c = 'PHYIMPL'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  引数の正当性のチェック
    !  Validate arguments
    !-----------------------------------------------------------------
    if (imax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'imax'
      goto 999
    end if
    if (jmax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'jmax'
      goto 999
    end if
    if (kmax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'kmax'
      goto 999
    end if
    if (DelTime < 0.0_DP) then
      stat = DC_ENEGATIVE
      cause_c = 'DelTime'
      goto 999
    end if

    
    !-----------------------------------------------------------------
    !  "PHYIMPL" の設定
    !  Configure the settings for "PHYIMPL"
    !-----------------------------------------------------------------
    phy_impl % imax    = imax   
    phy_impl % jmax    = jmax   
    phy_impl % kmax    = kmax   
    phy_impl % Grav    = Grav   
    phy_impl % Cp      = Cp     
    phy_impl % EL      = EL     
    phy_impl % DelTime = DelTime

    allocate( phy_impl % xy_SurfHeatCapacity (0:imax-1, 0:jmax-1) )
    allocate( phy_impl % xy_SurfCondition (0:imax-1, 0:jmax-1) )
    allocate( phy_impl % xy_GroundTempFlux (0:imax-1, 0:jmax-1) )
    phy_impl % xy_SurfHeatCapacity = xy_SurfHeatCapacity
    phy_impl % xy_SurfCondition = xy_SurfCondition
    phy_impl % xy_GroundTempFlux = xy_GroundTempFlux

    !-------------------------
    !  デフォルト値
    !  Default values
!!$    phy_impl % param_r = 0.0_DP
!!$    phy_impl % param_c = 'hogehoge'

    !-------------------------
    !  オプショナル引数からの値
    !  Values from optional arguments
!!$    phy_impl % param_i = param_i
!!$    if ( present(param_r) )  phy_impl % param_r = param_r
!!$    if ( present(param_c) )  phy_impl % param_c = param_c

    !-------------------------
    !  NAMELIST からの値
    !  Values from NAMELIST

!!$    if ( present_and_not_empty(nmlfile) ) then
!!$      call MessageNotify( 'M', subname, &
!!$        & 'Loading NAMELIST file "%c" ...', &
!!$        & c1 = trim(nmlfile) )
!!$      call NmlRead ( nmlfile = nmlfile, &      ! (in)
!!$        & param_i = phy_impl % param_i, &   ! (inout)
!!$        & param_r = phy_impl % param_r, &   ! (inout)
!!$        & param_c_ = phy_impl % param_c, &  ! (inout)
!!$        & err = err )                          ! (out)
!!$      if ( present_and_true(err) ) then
!!$        call MessageNotify( 'W', subname, &
!!$          & '"%c" can not be read.', &
!!$          & c1 = trim(nmlfile) )
!!$        stat = DC_ENOFILEREAD
!!$        cause_c = nmlfile
!!$        goto 999
!!$      end if
!!$    end if

    !-----------------------------------------------------------------
    !  設定値の正当性のチェック
    !  Validate setting values
    !-----------------------------------------------------------------
!!$    if ( phy_impl % param_i < 0 ) then
!!$      stat = DC_ENEGATIVE
!!$      cause_c = 'param_i'
!!$      goto 999
!!$    end if


    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
    phy_impl % initialized = .true.
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyImplicitCreate

  subroutine PhyImplicitGetMatrices( phy_impl, xyr_Press, xyr_UFlux, xyr_VFlux, xyr_TempFlux, xyr_QVapFlux, xyza_UVMatrix, xyza_TempMatrix, xyza_QVapMatrix, err )
    !
    ! 陰解行列の作成と取得を行います.
    !
    ! なお, 与えられた *phy_impl* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Create and get implicit matrices. 
    !
    ! If *phy_impl* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYIMPL), intent(inout):: phy_impl
    real(DP), intent(in):: xyr_Press (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! $ P_s $ . 地表面気圧 (半整数レベル). 
                              ! Surface pressure (half level)
    real(DP), intent(out):: xyr_UFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 東西風速フラックス. 
                              ! Zonal wind flux
    real(DP), intent(out):: xyr_VFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 南北風速フラックス. 
                              ! Meridional wind flux
    real(DP), intent(out):: xyr_TempFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 温度フラックス. 
                              ! Temperature flux
    real(DP), intent(out):: xyr_QVapFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(out):: xyza_UVMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1, -1:1)
                              ! 速度陰解行列. 
                              ! Implicit matrix about velocity 
    real(DP), intent(out):: xyza_TempMatrix  (0:phy_impl%imax-1, 0:phy_impl%jmax-1, -1:phy_impl%kmax-1, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP), intent(out):: xyza_QVapMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1, -1:1)
                              ! 比湿陰解行列. 
                              ! Implicit matrix about specific humidity
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer:: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level

    real(DP):: Grav      ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
    real(DP):: DelTime    ! $ \Delta t $ . タイムステップ. Time step

    real(DP):: xy_SurfHeatCapacity (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              ! 地表熱容量. 
                              ! Surface heat capacity
    integer:: xy_SurfCondition (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              ! 地表状態. 
                              ! Surface condition

    integer:: i, j, k         ! DO ループ用作業変数
                              ! Work variables for DO loop

    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyImplicitGetMatrices'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_impl % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYIMPL'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_impl* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_impl*
    !-----------------------------------------------------------------
    imax    = phy_impl % imax   
    jmax    = phy_impl % jmax   
    kmax    = phy_impl % kmax   
    Grav    = phy_impl % Grav   
    Cp      = phy_impl % Cp     
    DelTime = phy_impl % DelTime

    xy_SurfHeatCapacity = phy_impl % xy_SurfHeatCapacity
    xy_SurfCondition    = phy_impl % xy_SurfCondition

    !----------------------------------------------------------------
    !  陰解行列作成
    !  Create implicit matrices
    !----------------------------------------------------------------

    !-----------------------------------
    !  質量, 熱容量の項
    !  Mass and heat capacity terms
    xyza_UVMatrix  = 0.0_DP
    xyza_TempMatrix = 0.0_DP
    xyza_QVapMatrix = 0.0_DP

    do k = 0, kmax-1
      xyza_UVMatrix(:,:,k,0)  = ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Grav / ( 2.0_DP * DelTime )
      xyza_TempMatrix(:,:,k,0) = xyza_UVMatrix(:,:,k,0) * Cp
      xyza_QVapMatrix(:,:,k,0) = xyza_UVMatrix(:,:,k,0) * Cp
    end do

    do j = 0, jmax-1
      do i = 0, imax-1
        if ( xy_SurfCondition(i,j) >= 1 ) then 
          xyza_TempMatrix(i,j,-1,0) = xy_SurfHeatCapacity(i,j) / ( 2.0_DP * DelTime )
        else
          xyza_TempMatrix(i,j,-1,0) = 1.0_DP
        end if
      end do
    end do

    !-----------------------------------
    !  フラックスをリセット
    !  Reset fluxes
    xyr_UFlux    = 0.0_DP
    xyr_VFlux    = 0.0_DP
    xyr_TempFlux = 0.0_DP
    xyr_QVapFlux = 0.0_DP

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyImplicitGetMatrices

  logical function PhyImplicitInitialized( phy_impl ) result(result)
    !
    ! *phy_impl* が初期設定されている場合には .true. が,
    ! 初期設定されていない場合には .false. が返ります.
    !
    ! If *phy_impl* is initialized, .true. is returned.
    ! If *phy_impl* is not initialized, .false. is returned.
    !
    implicit none
    type(PHYIMPL), intent(in):: phy_impl
  continue
    result = phy_impl % initialized
  end function PhyImplicitInitialized

  subroutine PhyImplicitIntegrate( phy_impl, xyr_UFlux, xyr_VFlux, xyr_TempFlux, xy_SurfRadSFlux, xy_SurfRadLFlux, xyr_QVapFlux, xyza_UVMatrix, xyza_TempMatrix, xyza_QVapMatrix, xy_SurfUVMatrix, xyaa_SurfTempMatrix, xyaa_SurfQVapMatrix, xya_SurfRadLMatrix, xyz_DVerdiffUDt, xyz_DVerdiffVDt, xyz_DVerdiffTempDt, xy_DVerdiffSurfTempDt, xyz_DVerdiffQVapDt, err )
    !
    ! 陰解行列の作成と取得を行います.
    !
    ! なお, 与えられた *phy_impl* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Create and get implicit matrices. 
    !
    ! If *phy_impl* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYIMPL), intent(inout):: phy_impl
    real(DP), intent(in):: xyr_UFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 東西風速フラックス. 
                              ! Zonal wind flux
    real(DP), intent(in):: xyr_VFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 南北風速フラックス. 
                              ! Meridional wind flux
    real(DP), intent(in):: xyr_TempFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 温度フラックス. 
                              ! Temperature flux
    real(DP), intent(in):: xy_SurfRadSFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              ! 短波 (日射) フラックス (地表面). 
                              ! Short wave (insolation) flux on surface
    real(DP), intent(in):: xy_SurfRadLFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              ! 長波フラックス (地表面). 
                              ! Long wave flux on surface
    real(DP), intent(in):: xyr_QVapFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(in):: xyza_UVMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1, -1:1)
                              ! 速度陰解行列. 
                              ! Implicit matrix about velocity 
    real(DP), intent(in):: xyza_TempMatrix  (0:phy_impl%imax-1, 0:phy_impl%jmax-1, -1:phy_impl%kmax-1, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP), intent(in):: xyza_QVapMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1, -1:1)
                              ! 比湿陰解行列. 
                              ! Implicit matrix about specific humidity
    real(DP), intent(in):: xy_SurfUVMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              !  速度陰解行列: 地表. 
                              ! Implicit matrix about velocity: surface
    real(DP), intent(in):: xyaa_SurfTempMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:1, -1:1)

                              ! 温度陰解行列: 地表. 
                              ! Implicit matrix about temperature: surface
    real(DP), intent(in):: xyaa_SurfQVapMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:1, -1:1)
                              ! 比湿陰解行列: 地表. 
                              ! Implicit matrix about specific humidity: surface
    real(DP), intent(in):: xya_SurfRadLMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, -1:1)
                              ! $ T $ 陰解行列: 地表. 
                              ! $ T $ implicit matrix: surface
    real(DP), intent(out):: xyz_DVerdiffUDt (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1)
                              ! 東西成分 鉛直拡散加速度. 
                              ! Zonal acceleration by vertical diffusion
    real(DP), intent(out):: xyz_DVerdiffVDt (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1)
                              ! 南北成分 鉛直拡散加速度. 
                              ! Meridional acceleration by vertical diffusion
    real(DP), intent(out):: xyz_DVerdiffTempDt (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1)
                              ! 鉛直拡散加熱率. 
                              ! Temperature tendency by vertical diffusion
    real(DP), intent(out):: xy_DVerdiffSurfTempDt (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              ! 地表面 鉛直拡散加熱率. 
                              ! Surface temperature tendency by vertical diffusion
    real(DP), intent(out):: xyz_DVerdiffQVapDt (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1)
                              ! 鉛直拡散加湿率. 
                              ! Specific humidity tendency by vertical diffusion
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer:: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level

    real(DP):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
    real(DP):: EL        ! $ L $ .      水の凝結の潜熱. Latent heat of condensation of water vapor
    real(DP):: DelTime    ! $ \Delta t $ . タイムステップ. Time step

    integer:: xy_SurfCondition (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              ! 地表状態. 
                              ! Surface condition
    real(DP):: xy_GroundTempFlux (0:phy_impl%imax-1, 0:phy_impl%jmax-1)
                              ! 地中熱フラックス. 
                              ! Ground temperature flux

    real(DP):: xyz_DelTempQVap (0:phy_impl%imax-1, 0:phy_impl%jmax-1, -phy_impl%kmax:phy_impl%kmax)
                              ! $ T q $ の時間変化. 
                              ! Tendency of $ T q $ 
    real(DP):: xyza_TempQVapLUMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, -phy_impl%kmax:phy_impl%kmax, -1:1)
                              ! LU 行列. 
                              ! LU matrix
    real(DP):: xyza_UVLUMatrix (0:phy_impl%imax-1, 0:phy_impl%jmax-1, 0:phy_impl%kmax-1,-1:1)
                              ! LU 行列. 
                              ! LU matrix

    integer:: i, j, k, l      ! DO ループ用作業変数
                              ! Work variables for DO loop

    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyImplicitIntegrate'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_impl % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYIMPL'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_impl* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_impl*
    !-----------------------------------------------------------------
    imax    = phy_impl % imax   
    jmax    = phy_impl % jmax   
    kmax    = phy_impl % kmax   
    Cp      = phy_impl % Cp     
    EL      = phy_impl % EL     
    DelTime = phy_impl % DelTime

    xy_SurfCondition    = phy_impl % xy_SurfCondition
    xy_GroundTempFlux   = phy_impl % xy_GroundTempFlux

    !----------------------------------------------------------------
    !  東西風速, 南北風速の計算
    !  Calculate zonal and meridional wind
    !----------------------------------------------------------------
    xyza_UVLUMatrix = xyza_UVMatrix
    xyza_UVLUMatrix(:,:,0,0)  = xyza_UVLUMatrix(:,:,0,0) + xy_SurfUVMatrix
    
    call PhyImplicitLUDecomp3( xyza_UVLUMatrix, imax * jmax, kmax ) ! (in)
    
    do k = 0, kmax - 1
      xyz_DVerdiffUDt(:,:,k) = xyr_UFlux(:,:,k) - xyr_UFlux(:,:,k+1)
      xyz_DVerdiffVDt(:,:,k) = xyr_VFlux(:,:,k) - xyr_VFlux(:,:,k+1)
    end do

    call PhyImplicitLUSolve3( xyz_DVerdiffUDt, xyza_UVLUMatrix, 1, imax * jmax, kmax ) ! (in)
    
    call PhyImplicitLUSolve3( xyz_DVerdiffVDt, xyza_UVLUMatrix, 1, imax * jmax, kmax ) ! (in)
    
    !----------------------------------------------------------------
    !  温度と比湿の計算
    !  Calculate temperature and specific humidity
    !----------------------------------------------------------------
    do l = -1, 1
      
      do k = 1, kmax
        xyza_TempQVapLUMatrix(:,:,k,l)   = xyza_TempMatrix(:,:,k-1,l)
        xyza_TempQVapLUMatrix(:,:,-k,-l) = xyza_QVapMatrix(:,:,k-1,l)
      end do
      
      xyza_TempQVapLUMatrix(:,:,1,l) = xyza_TempMatrix(:,:,0,l) + xyaa_SurfTempMatrix(:,:,1,l)

      xyza_TempQVapLUMatrix(:,:,-1,-l) = xyza_QVapMatrix(:,:,0,l) + xyaa_SurfQVapMatrix(:,:,1,l)
      
    end do
    
    xyza_TempQVapLUMatrix(:,:,0,0) = xyza_TempMatrix(:,:,-1,0) + xyaa_SurfTempMatrix(:,:,0,0) + xyaa_SurfQVapMatrix(:,:,0,0) + xya_SurfRadLMatrix(:,:,0)
    
    xyza_TempQVapLUMatrix(:,:,0,1) = xyaa_SurfTempMatrix(:,:,0,1) + xya_SurfRadLMatrix(:,:,1)
    
    xyza_TempQVapLUMatrix(:,:,0,-1) =  xyaa_SurfQVapMatrix(:,:,0,1) 
    
    call PhyImplicitLUDecomp3( xyza_TempQVapLUMatrix, imax * jmax, (2 * kmax) + 1 ) ! (in)
    
    do k = 1, kmax
      xyz_DelTempQVap(:,:,k)  = xyr_TempFlux(:,:,k-1) - xyr_TempFlux(:,:,k)
      xyz_DelTempQVap(:,:,-k) = xyr_QVapFlux(:,:,k-1) - xyr_QVapFlux(:,:,k)
    end do
    
    xyz_DelTempQVap(:,:,0) = - xy_SurfRadSFlux - xy_SurfRadLFlux - xyr_TempFlux(:,:,0) - xyr_QVapFlux(:,:,0) + xy_GroundTempFlux
    
    call PhyImplicitLUSolve3( xyz_DelTempQVap, xyza_TempQVapLUMatrix, 1, imax * jmax , (2 * kmax) + 1 ) ! (in)
    
    !----------------------------------------------------------------
    !  時間変化率の計算
    !  Calculate tendency
    !----------------------------------------------------------------
    do k = 1, kmax
      xyz_DVerdiffUDt(:,:,k-1) = xyz_DVerdiffUDt(:,:,k-1) / ( 2.0_DP * DelTime )
      xyz_DVerdiffVDt(:,:,k-1) = xyz_DVerdiffVDt(:,:,k-1) / ( 2.0_DP * DelTime )
      xyz_DVerdiffTempDt(:,:,k-1) = xyz_DelTempQVap(:,:,k) / ( 2.0_DP * DelTime )
      xyz_DVerdiffQVapDt(:,:,k-1) = xyz_DelTempQVap(:,:,-k) / ( 2.0_DP * DelTime ) / EL * Cp
    end do

    do j = 0, jmax-1
      do i = 0, imax-1
        if ( xy_SurfCondition(i,j) >= 1 ) then
          xy_DVerdiffSurfTempDt(i,j) = xyz_DelTempQVap(i,j,0) / ( 2.0_DP * DelTime )
        else
          xy_DVerdiffSurfTempDt(i,j) = 0.0_DP
        end if
      end do
    end do

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyImplicitIntegrate

  subroutine PhyImplicitLUDecomp3( jna_LUMatrix, JDimMax, NDimMax )
    !
    ! 3 重対角行列の LU 分解を行います. 
    !
    ! なお, 与えられた *phy_impl* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! LU decomposition of triple diagonal matrix.
    !
    ! If *phy_impl* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    integer, intent(in):: JDimMax
    integer, intent(in):: NDimMax
    real(DP), intent(inout):: jna_LUMatrix(JDimMax, NDimMax, -1:1)
                              ! LU 行列. 
                              ! LU matrix

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: j, n            ! DO ループ用作業変数
                              ! Work variables for DO loop

    character(*), parameter:: subname = 'PhyImplicitLUDecomp3'
  continue
    call BeginSub( subname )

    !----------------------------------------------------------------
    !  LU 分解
    !  LU decomposition
    !----------------------------------------------------------------
    do j = 1, JDimMax
      jna_LUMatrix(j,1,1) = jna_LUMatrix(j,1,1) / jna_LUMatrix(j,1,0)
    end do

    do n = 2, NDimMax-1
      do j = 1, JDimMax
        jna_LUMatrix(j,n,0) = jna_LUMatrix(j,n,0) - jna_LUMatrix(j,n,-1) * jna_LUMatrix(j,n-1,1)

        jna_LUMatrix(j,n,1) = jna_LUMatrix(j,n,1) /jna_LUMatrix(j,n,0)
      end do
    end do

    do j = 1, JDimMax
      jna_LUMatrix(j,NDimMax,0) = jna_LUMatrix(j,NDimMax,0) - jna_LUMatrix(j,NDimMax,-1) * jna_LUMatrix(j,NDimMax-1,1)
    end do

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
    call EndSub( subname )
  end subroutine PhyImplicitLUDecomp3

  subroutine PhyImplicitLUSolve3( ijn_Vector, jna_LUMatrix, IDimMax, JDimMax, NDimMax )
    !
    ! LU 分解による解の計算 (3重対角行列用)
    !
    ! Solve with LU decomposition (For triple diagonal matrix). 
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    integer, intent(in):: IDimMax
    integer, intent(in):: JDimMax
    integer, intent(in):: NDimMax
    real(DP), intent(in):: jna_LUMatrix(JDimMax, NDimMax, -1:1)
                              ! LU 行列. 
                              ! LU matrix
    real(DP), intent(inout):: ijn_Vector(IDimMax, JDimMax, NDimMax)
                              ! 右辺ベクトル / 解
                              ! Right-hand side vector / solution

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: i, j, n         ! DO ループ用作業変数
                              ! Work variables for DO loop

    character(*), parameter:: subname = 'PhyImplicitLUSolve3'
  continue
    call BeginSub( subname )

    !-----------------------------------------------------------------
    !  前進代入
    !  Forward substitution
    !-----------------------------------------------------------------
    do i = 1, IDimMax
      do j = 1, JDimMax
        ijn_Vector(i,j,1) = ijn_Vector(i,j,1) / jna_LUMatrix(j,1,0)
      end do
    end do

    do n = 2, NDimMax
      do i = 1, IDimMax
        do j = 1, JDimMax
          ijn_Vector(i,j,n) = (   ijn_Vector(i,j,n) - ijn_Vector(i,j,n-1) * jna_LUMatrix(j,n,-1) ) / jna_LUMatrix(j,n,0)
        end do
      end do
    end do

    !-----------------------------------------------------------------
    !  後退代入
    !  Backward substitution
    !-----------------------------------------------------------------
    do n = NDimMax-1, 1, -1
      do i = 1, IDimMax
        do j = 1, JDimMax
          ijn_Vector(i,j,n) = ijn_Vector(i,j,n) - ijn_Vector(i,j,n+1) * jna_LUMatrix(j,n,1)
        end do
      end do
    end do

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
    call EndSub( subname )
  end subroutine PhyImplicitLUSolve3

  subroutine PhyImplicitNmlRead( nmlfile, err )
    !
    ! NAMELIST ファイル *nmlfile* から値を入力するための
    ! 内部サブルーチンです. Create 内で呼び出されることを
    ! 想定しています.
    !
    ! 値が NAMELIST ファイル内で指定されていない場合には,
    ! 入力された値がそのまま返ります.
    !
    ! なお, *nmlfile* に空文字が与えられた場合, または
    ! 与えられた *nmlfile* を読み込むことができない場合, 
    ! プログラムはエラーを発生させます.
    !
    ! This is an internal subroutine to input values from 
    ! NAMELIST file *nmlfile*. This subroutine is expected to be
    ! called by "Create".
    !
    ! A value not specified in NAMELIST file is returned
    ! without change.
    !
    ! If *nmlfile* is empty, or *nmlfile* can not be read, 
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_iounit, only: FileOpen
    use dc_message, only: MessageNotify
    use dc_present, only: present_and_true
    use dc_error, only: StoreError, DC_NOERR, DC_ENOFILEREAD
    implicit none
    character(*), intent(in):: nmlfile
                              ! NAMELIST ファイルの名称. 
                              ! NAMELIST file name
!!$    integer, intent(inout):: param_i
!!$    real(DP), intent(inout):: param_r
!!$    character(*), intent(inout):: param_c_
!!$    character(TOKEN):: param_c
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

!!$    namelist /phy_implicit_nml/ &
!!$      & param_i, param_r, param_c
                              ! phy_implicit モジュール用
                              ! NAMELIST 変数群名.
                              !
                              ! phy_implicit#Create を使用する際に, 
                              ! オプショナル引数 *nmlfile* へ NAMELIST 
                              ! ファイル名を指定することで, そのファイルから
                              ! この NAMELIST 変数群を読み込みます.
                              !
                              ! NAMELIST group name for 
                              ! "phy_implicit" module.
                              ! 
                              ! If a NAMELIST filename is specified to 
                              ! an optional argument *nmlfile* 
                              ! when "phy_implicit#Create" is used, 
                              ! this NAMELIST group is loaded from 
                              ! the file.

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
!!$    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
!!$                              ! IOSTAT of NAMELIST read
    character(*), parameter:: subname = 'PhyImplicitNmlRead'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''



    !-----------------------------------------------------------------
    !  文字型引数を NAMELIST 変数群へ代入
    !  Substitute character arguments to NAMELIST group
    !-----------------------------------------------------------------
!!$    param_c = param_c_

    !----------------------------------------------------------------
    !  NAMELIST ファイルのオープン
    !  Open NAMELIST file
    !----------------------------------------------------------------
    call FileOpen( unit = unit_nml, file = nmlfile, mode = 'r', err = err )                   ! (out)
    if ( present_and_true(err) ) then
      stat = DC_ENOFILEREAD
      cause_c = nmlfile
      goto 999
    end if


    !-----------------------------------------------------------------
    !  NAMELIST 変数群の取得
    !  Get NAMELIST group
    !-----------------------------------------------------------------
!!$    read( unit = unit_nml, & ! (in)
!!$      & nml = phy_implicit_nml, iostat = iostat_nml ) ! (out)
!!$    if ( iostat_nml == 0 ) then
!!$      call MessageNotify( 'M', subname, &
!!$        & 'NAMELIST group "%c" is loaded from "%c".', &
!!$        & c1 = 'phy_implicit_nml', c2 = trim(nmlfile) )
!!$      write(STDOUT, nml = phy_implicit_nml)
!!$    else
!!$      call MessageNotify( 'W', subname, &
!!$        & 'NAMELIST group "%c" is not found in "%c" (iostat=%d).', &
!!$        & c1 = 'phy_implicit_nml', c2 = trim(nmlfile), &
!!$        & i = (/iostat_nml/) )
!!$    end if

    close( unit_nml )

    !-----------------------------------------------------------------
    !  NAMELIST 変数群を文字型引数へ代入
    !  Substitute NAMELIST group to character arguments
    !-----------------------------------------------------------------
!!$    param_c_ = param_c

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyImplicitNmlRead

  subroutine PhyImplicitPutLine( phy_impl, unit, indent, err )
    !
    ! 引数 *phy_impl* に設定されている情報を印字します.
    ! デフォルトではメッセージは標準出力に出力されます. 
    ! *unit* に装置番号を指定することで, 出力先を変更することが可能です.
    !
    ! Print information of *phy_impl*.
    ! By default messages are output to standard output.
    ! Unit number for output can be changed by *unit* argument.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYIMPL), intent(in):: phy_impl
    integer, intent(in), optional:: unit
                              ! 出力先の装置番号.
                              ! デフォルトの出力先は標準出力.
                              !
                              ! Unit number for output.
                              ! Default value is standard output.
    character(*), intent(in), optional:: indent
                              ! 表示されるメッセージの字下げ.
                              !
                              ! Indent of displayed messages.
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    integer:: out_unit
    integer:: indent_len
    character(STRING):: indent_str
    character(*), parameter:: subname = 'PhyImplicitPutLine'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( present(unit) ) then
      out_unit = unit
    else
      out_unit = STDOUT
    end if

    indent_len = 0
    indent_str = ''
    if ( present(indent) ) then
      if ( len(indent) /= 0 ) then
        indent_len = len(indent)
        indent_str(1:indent_len) = indent
      end if
    end if


    !-----------------------------------------------------------------
    !  "PHYIMPL" の設定の印字
    !  Print the settings for "PHYIMPL"
    !-----------------------------------------------------------------
    if ( phy_impl % initialized ) then
      call Printf( out_unit, indent_str(1:indent_len) // '#<PHYIMPL:: @initialized=%y', l = (/phy_impl % initialized/) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @imax=%d @jmax=%d @kmax=%d', i = (/phy_impl % imax, phy_impl % jmax, phy_impl % kmax/) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @Grav=%f @Cp=%f @EL=%f', d = (/ phy_impl % Grav, phy_impl % Cp, phy_impl % EL /) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @DelTime=%f', d = (/ phy_impl % DelTime /) )

      call PutLine( phy_impl % xy_SurfHeatCapacity, unit = out_unit, lbounds = lbound(phy_impl % xy_SurfHeatCapacity), ubounds = ubound(phy_impl % xy_SurfHeatCapacity), indent = indent_str(1:indent_len) // ' @xy_SurfHeatCapacity=' )

      call PutLine( phy_impl % xy_SurfCondition, unit = out_unit, lbounds = lbound(phy_impl % xy_SurfCondition), ubounds = ubound(phy_impl % xy_SurfCondition), indent = indent_str(1:indent_len) // ' @xy_SurfCondition=' )

      call PutLine( phy_impl % xy_GroundTempFlux, unit = out_unit, lbounds = lbound(phy_impl % xy_GroundTempFlux), ubounds = ubound(phy_impl % xy_GroundTempFlux), indent = indent_str(1:indent_len) // ' @xy_GroundTempFlux=' )

      call Printf( out_unit, indent_str(1:indent_len) // '>' )
    else
      call Printf( out_unit, indent_str(1:indent_len) // '#<PHYIMPL:: @initialized=%y>', l = (/phy_impl % initialized/) )
    end if

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyImplicitPutLine