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

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)   :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速
    real(DP), intent(in)   :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速
    real(DP), intent(in)   :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
                                           !温位の擾乱成分    
    real(DP), intent(in)   :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
                                           !温位の擾乱成分    
    real(DP), intent(in)   :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                           !温位の擾乱成分    
    real(DP), intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)               :: pyz_DVelXDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xqz_DVelYDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)               :: pyz_Momflux(imin:imax,jmin:jmax,kmin:kmax)
                                           !運動量フラックス
    real(DP)               :: xqz_Momflux(imin:imax,jmin:jmax,kmin:kmax)
                                           !運動量フラックス
    real(DP)               :: xyz_Heatflux(imin:imax,jmin:jmax,kmin:kmax)
                                           !地表面熱フラックス
    real(DP)               :: xyzf_QMixflux(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                           !物質的フラックス
    real(DP)               :: pyz_VelXSfc(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xqz_VelYSfc(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_VelXSfc(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: xyz_VelYSfc(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: TempSfc 
    real(DP)               :: PressSfc 
    real(DP)               :: xyz_VelSfc(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    integer                :: kz            !配列添字
    integer                :: s             !ループ変数

    ! 初期化
    !
    kz = 1
    xyz_Heatflux = 0.0d0
    xyzf_QMixflux = 0.0d0
    pyz_Momflux = 0.0d0
    xqz_Momflux = 0.0d0

    pyz_DVelXDt0  = pyz_DVelXDt
    xqz_DVelYDt0  = xqz_DVelYDt
    xyz_DPTempDt0 = xyz_DPTempDt
    xyzf_DQMixDt0 = xyzf_DQMixDt

    xyz_TempAll  = (xyz_Exner + xyz_ExnerBZ) * (xyz_PTemp + xyz_PTempBZ)
    xyzf_QMixAll = xyzf_QMix + xyzf_QMixBZ
    
    pyz_VelXSfc = pyz_VelX + VelX0
    xqz_VelYSfc = xqz_VelY + VelY0    

    ! 地表面での値
    ! 
    TempSfc = xyz_TempBZ(1,1,kz)  !水平一様なので, i=j=1 で代表させる
    PressSfc= xyz_PressBZ(1,1,kz) !水平一様なので, i=j=1 で代表させる
    
    !地表面熱フラックスによる加熱率を計算
    !  * 単位は K/s
    !  * エクスナー関数は基本場の値で代表させる.     
    !  * 格子点 xz では, 物理領域の最下端の添え字は kz = 1
    xyz_VelXSfc = xyz_avr_pyz(pyz_VelXSfc)
    xyz_VelYSfc = xyz_avr_xqz(xqz_VelYSfc)    
    
    xyz_VelSfc = SQRT(  xyz_VelXSfc ** 2.0 + xyz_VelYSfc ** 2.0 )

    xyz_Heatflux(:,:,kz) = MAX( 0.0d0, - Bulk * xyz_VelSfc(:,:,kz) * ( xyz_TempAll(:,:,kz) - TempSfc ) / ( z_dz(kz) * 5.0d-1 ) )

    do s = 1, CondNum
      xyzf_QMixflux(:,:,kz,IdxCG(s)) = MAX( 0.0d0, - Bulk * xyz_VelSfc(:,:,kz) * (  xyzf_QMixAll(:,:,kz,s) - SvapPress( SpcWetID(IdxCC(s)), TempSfc ) / PressSfc * (MolWtWet(IdxCG(s)) / MolWtDry) ) / ( z_dz(kz) * 5.0d-1 ) )
    end do

    !地表面運動量フラックスによる変化率を計算
    !  * 単位は m/s^2
    !  * 格子点 xz では, 物理領域の最下端の添え字は kz = 1
    pyz_MomFlux(:,:,kz) = MAX( 0.0d0, - Bulk * abs(pyz_VelXSfc(:,:,kz)) * pyz_VelX(:,:,kz) / ( z_dz(kz) * 5.0d-1 ) )

    xqz_MomFlux(:,:,kz) = MAX( 0.0d0, - Bulk * abs(xyz_VelYSfc(:,:,kz)) * xqz_VelY(:,:,kz) / ( z_dz(kz) * 5.0d-1 ) )

    xyz_DPTempDt = xyz_DPTempDt0 + xyz_Heatflux
    xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Qmixflux
    pyz_DVelXDt  = pyz_DVelXDt0  + pyz_MomFlux
    xqz_DVelYDt  = xqz_DVelYDt0  + xqz_MomFlux

    call HistoryAutoPut(TimeN, 'PTempSfc', xyz_HeatFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSfc',  pyz_MomFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSfc',  xqz_MomFlux(1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Sfc', xyzf_Qmixflux(1:nx,1:ny,1:nz,s))
    end do
  end subroutine Surfaceflux_Bulk_forcing

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

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer    :: l, unit

    !---------------------------------------------------------------    
    ! NAMELIST から情報を取得
    !
    NAMELIST /surfaceflux_bulk_nml/ Bulk, VelX0, VelY0
    
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=surfaceflux_bulk_nml)
    close(unit)  

    if (myrank == 0) then 
      call MessageNotify( "M", "SurfaceSfc", "Bulk = %f", d=(/Bulk/))
      call MessageNotify( "M", "SurfaceSfc", "VelX0 = %f", d=(/VelX0/))
      call MessageNotify( "M", "SurfaceSfc", "VelY0 = %f", d=(/VelY0/))
    end if

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

    call HistoryAutoAddVariable( varname='VelXSfc', dims=(/'x','y','z','t'/), longname='flux of heating from surface', units='kg.m.s-1', xtype='double')

    call HistoryAutoAddVariable( varname='VelYSfc', dims=(/'x','y','z','t'/), longname='flux of momentum from surface', units='kg.m.s-1', xtype='double')

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

  end subroutine Surfaceflux_Bulk_init