function pz_MomFluxBulk( Vel )
    !
    ! 下部境界からのフラックスによる運動量の変化率を,
    ! バルク方法に基づいて計算する.
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in)   :: Vel(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !水平風速
    real(8)               :: pz_MomFluxBulk(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !地表面熱フラックス
    integer               :: kz            !配列添字


    !初期化
    !  * 全ての値をゼロに固定
    pz_MomFluxBulk = 0.0d0

    !地表面運動量フラックスによる変化率を計算
    !  * 単位は m/s^2
    !  * 格子点 xz では, 物理領域の最下端の添え字は RegZMin+1

    kz = RegZMin+1

    pz_MomFluxBulk(:,kz) = - Bulk * abs(Vel(:,kz)) * Vel(:,kz) /DelZ

  end function pz_MomFluxBulk

  function xz_HeatFluxBulk( xz_PotTemp, pz_VelX )
    ! 
    ! 下部境界からのフラックスによる温度の変化率を,
    ! バルク方法に基づいて計算する.
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)   :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !温位の擾乱成分    
    real(8), intent(in)   :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !水平風速
    real(8)               :: xz_HeatFluxBulk(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !地表面熱フラックス
    real(8)               :: VelX = 3.0d0  !下層での水平速度嵩上げ値
    integer               :: kz            !配列添字

    real(8), allocatable  :: xz_VelX(:,:)  !水平風速 (xz 格子)


    !作業配列の割り付け
    allocate(xz_VelX(DimXMin:DimXMax,DimZMin:DimZMax))


    !初期化
    !  * 全ての値をゼロに固定
    xz_HeatFluxBulk = 0.0d0
    
    !地表面熱フラックスによる加熱率を計算
    !  * 単位は K/s
    !  * エクスナー関数は基本場の値で代表させる.     
    !  * 格子点 xz では, 物理領域の最下端の添え字は RegZMin+1

    kz = RegZMin+1
    xz_VelX = xz_avr_pz(pz_VelX)

    xz_HeatFluxBulk(:,kz) = MAX( 0.0d0, - Bulk * SQRT(xz_VelX(:,kz)**2.0 + VelX**2.0) * (xz_ExnerBasicZ(:,kz) * (xz_PotTemp(:,kz)  + xz_PotTempBasicZ(:,kz)) - TempSfc )/DelZ )

    !作業配列の解放
    deallocate(xz_VelX)
    
  end function xz_HeatFluxBulk

  function xz_MixRtFluxBulk( xz_MixRt, pz_VelX, xz_Exner, ID )
    !
    ! 下部境界からのフラックスによる凝結成分混合比の変化率を,
    ! バルク方法に基づいて計算する.
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)   :: xz_MixRt(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !温位の擾乱成分    
    real(8), intent(in)   :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !水平風速
    real(8), intent(in)   :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !圧力の擾乱成分    
    integer, intent(in)   :: ID            !凝結成分の ID 番号
    real(8)               :: xz_MixRtFluxBulk(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !地表面混合比フラックス
    real(8)               :: VelX = 3.0d0  !下層での水平速度嵩上げ値
    integer               :: kz            !配列添字

    real(8), allocatable  :: xz_VelX(:,:)  !水平風速 (xz 格子)


    !作業配列の割り付け
    allocate(xz_VelX(DimXMin:DimXMax,DimZMin:DimZMax))


    !初期化
    !  * 全ての値をゼロに固定
    xz_MixRtFluxBulk = 0.0d0

    kz = RegZMin+1
    xz_VelX = xz_avr_pz(pz_VelX)
    
    xz_MixRtFluxBulk(:,kz) = max( 0.0d0, - Bulk * SQRT(xz_VelX(:,kz)**2.0 + VelX**2.0) * (xz_MixRt(:,kz) - SvapPress( ID, TempSfc ) / ( ( xz_Exner(:,kz) + xz_ExnerBasicZ(:,kz) ) **( CpDry / GasRDry )*PressSfc ) )/DelZ )

    !作業配列の解放
    deallocate(xz_VelX)

    
  end function xz_MixRtFluxBulk