subroutine Surfaceflux_Diff_forcing( xyz_PTemp, xyzf_QMix, xyz_DPTempDt, xyz_DExnerDt, xyzf_DQMixDt )
    ! 
    ! 下部境界からのフラックスによる温度の変化率を,
    ! バルク方法に基づいて計算する.
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)   :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
                                           !温位の擾乱成分    
    real(DP), intent(in)   :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                           !温位の擾乱成分    
    real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)               :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_DExnerDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)               :: xyz_Heatflux(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_Exnerflux(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyzf_QMixflux(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                            !地表面熱フラックス
    integer, parameter     :: kz = 1        !配列添字
    integer                :: l             !ループ変数

    ! 初期化
    !
    xyz_HeatFlux  = 0.0d0
    xyz_ExnerFlux = 0.0d0
    xyzf_QMixFlux = 0.0d0

    xyz_DPTempDt0 = xyz_DPTempDt
    xyz_DExnerDt0 = xyz_DExnerDt
    xyzf_DQMixDt0 = xyzf_DQMixDt

    !地表面熱フラックスによる加熱率を計算
    !  * 単位は K/s
    !  * エクスナー関数は基本場の値で代表させる.     
    !  * 格子点 xz では, 物理領域の最下端の添え字は kz = 1

    xyz_HeatFlux(:,:,kz) = - Kappa * xyz_PTemp(:,:,kz) * xyz_ExnerBZ(:,:,kz) / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) 

    do l = 1, GasNum
      xyzf_QMixFlux(:,:,kz,l) = max( 0.0d0, - Kappa * xyzf_QMix(:,:,kz,l) / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) )
    end do

    xyz_DPTempDt = xyz_DPTempDt0 + xyz_Heatflux
    xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Qmixflux

    if ( .not. FlagDExnerDtSurf ) then
      xyz_ExnerFlux = xyz_DExnerDt_xyz_xyzf( xyz_HeatFlux, xyzf_QMixFlux )
    else 
      xyz_ExnerFlux = 0.0d0
    end if
    xyz_DExnerDt = xyz_DExnerDt0 + xyz_ExnerFlux

    call HistoryAutoPut(TimeN, 'PTempFlux', xyz_HeatFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerFlux', xyz_ExnerFlux(1:nx,1:ny,1:nz))
    do l = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Flux', xyzf_Qmixflux(1:nx,1:ny,1:nz,l))
    end do    

    ! Set margin
    !
!    call SetMargin_xyz(xyz_DPTempDt)
!    call SetMargin_xyzf(xyzf_DQMixDt)

  end subroutine Surfaceflux_Diff_forcing

  subroutine Surfaceflux_Diff_init
    !
    !NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 
    !

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer    :: l, unit

    !---------------------------------------------------------------    
    ! NAMELIST から情報を取得
    !
    NAMELIST /surfaceflux_diff_nml/ Kappa, FlagDExnerDtSurf 

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=surfaceflux_diff_nml)
    close(unit)  

    if (myrank == 0) then 
      call MessageNotify( "M", module_name, "Kappa = %f", d=(/Kappa/))
      call MessageNotify( 'M', module_name, "FlagDExnerDtSurf = %b", l = (/ FlagDExnerDtSurf /) )
    end if

    call HistoryAutoAddVariable( varname='PTempFlux', dims=(/'x','y','z','t'/), longname='surface flux of potential temperature', units='kg.kg-1.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerFlux', dims=(/'x','y','z','t'/), longname='surface flux of Exner function', units='kg.kg-1.s-1', xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_Flux', dims=(/'x','y','z','t'/), longname='Surface Flux term of ' //trim(SpcWetSymbol(l))//' mixing ratio', units='kg.kg-1.s-1', xtype='float')
    end do

  end subroutine Surfaceflux_Diff_init