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

    ! 暗黙の型宣言禁止
    ! Implicit none
    !
    implicit none

    ! 変数
    ! variables
    !
    real(DP), intent(in)   :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速
                              ! X-component velocity
    real(DP), intent(in)   :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速
                              ! Y-component velocity
    real(DP), intent(in)   :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温位
                              ! Potential temperature
    real(DP), intent(in)   :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力関数
                              ! Exner function
    real(DP), intent(in)   :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 混合比
                              ! Mixing ration
    real(DP), intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 風速時間変化率
                              ! X-component velocity tendency
    real(DP), intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 風速時間変化率
                              ! Y-component velocity tendency
    real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温位時間変化率
                              ! Potential tempreture tendency
    real(DP), intent(inout):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力関数時間変化率
                              ! Exner function tendency
    real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 混合比時間変化率
                              ! Mixing ratio tendency

    ! 作業変数
    ! Work variables
    real(DP) :: xy_SurfBulkRiNum(imin:imax,jmin:jmax)
                              ! バルクリチャードソン数
                              ! Bulk Richardson number
    real(DP) :: xy_SurfRoughLength(imin:imax,jmin:jmax)
                              ! 祖度長さ
                              ! Roughness length
    real(DP) :: xy_SurfVelBulkCoef(imin:imax,jmin:jmax)
                              ! バルク係数(運動量)
                              ! Bulk coefficient for momentum
    real(DP) :: xy_SurfTempBulkCoef(imin:imax,jmin:jmax)
                              ! バルク係数(熱)
                              ! Bulk coefficient for heat
    real(DP) :: xy_SurfQmixBulkCoef(imin:imax,jmin:jmax)
                              ! バルク係数(混合比)
                              ! Bulk coefficient for mixing ratio
    real(DP) :: py_VelXflux (imin:imax,jmin:jmax)
                              ! x 方向速度フラックス
                              ! velocity flux in x
    real(DP) :: xq_VelYflux (imin:imax,jmin:jmax)
                              ! y 方向速度フラックス
                              ! celocity flux in y
    real(DP) :: xy_PTempFlux(imin:imax,jmin:jmax)
                              ! 温位フラックス
                              ! potential temperature flux
    real(DP) :: xy_ExnerFlux(imin:imax,jmin:jmax)
                              !
                              !
    real(DP) :: xyf_QMixFlux(imin:imax,jmin:jmax,ncmax)
                              ! 凝結成分混合比フラックス
                              ! Mixing ratio flux
    real(DP) :: xyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速 (xyz 格子)
                              ! X-component velocity (xyz grid)
    real(DP) :: xyz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速 (xyz 格子)
                              ! Y-component velocity (xyz grid)
    real(DP) :: xyz_AbsVel(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速絶対値 (xyz 格子)
                              ! Absolute value of horizontal velocity (xyz grid)
    real(DP) :: pyz_AbsVel(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速絶対値 (pyz 格子)
                              ! Absolute value of horizontal velocity (pyz grid)
    real(DP) :: xqz_AbsVel(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速絶対値 (xqz 格子)
                              ! Absolute value of horizontal velocity (xqz grid)
    real(DP) :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温位(基本場 + 擾乱)
                              ! Total value of potential temperature
    real(DP) :: xyz_ExnerAll (imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力関数(基本場 + 擾乱)
                              ! Total value of exner function
    real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 混合比(基本場 + 擾乱)
                              ! Total value of mixing ratios
    real(DP) :: xy_DPTempDtBulk(imin:imax,jmin:jmax)
                              ! 時間変化(温位)
                              ! potential temperature tendency by surface flux
    real(DP) :: xy_DExnerDtBulk(imin:imax,jmin:jmax)
                              ! 時間変化(圧力関数)
                              ! Exner function tendency by surface flux
    real(DP) :: xyf_DQMixDtBulk(imin:imax,jmin:jmax, ncmax)
                              ! 時間変化(混合比)
                              ! Mixing ratio tendency by surface flux
    real(DP) :: py_DVelXDtBulk (imin:imax,jmin:jmax)
                              ! 時間変化(U)
                              ! x-component velocity tendency by surface flux
    real(DP) :: xq_DVelYDtBulk (imin:imax,jmin:jmax)
                              ! 時間変化(V)
                              ! y-component velocity tendency by surface flux
    real(DP) :: xyz_DPTempDtBulk(imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(温位)
                              ! potential temperature tendency by surface flux
    real(DP) :: xyz_DExnerDtBulk(imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(圧力関数)
                              ! Exner function tendency by surface flux
    real(DP) :: xyzf_DQMixDtBulk(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 時間変化(混合比)
                              ! Mixing ratio tendency by surface flux
    real(DP) :: pyz_DVelXDtBulk (imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(U)
                              ! x-component velocity tendency by surface flux
    real(DP) :: xqz_DVelYDtBulk (imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(V)
                              ! y-component velocity tendency by surface flux
    real(DP) :: ExnerBZSfc    ! 地表面圧力関数 
                              ! Basic state Exner function at the surface
    real(DP) :: xy_PressSfc(imin:imax,jmin:jmax)
                              ! Total pressure at the surface
    integer  :: kz            ! 配列添字
                              ! Arrzy index
    integer  :: s             ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituen                         

    ! 初期化
    ! Initialization
    ! 
    kz = 1

    ! 祖度長さの指定
    ! Specify surface length
    xy_SurfRoughLength = SfcRoughLength

    ! 全量の計算
    ! Calculate total value of thermodynamic variables
    ! 
    xyz_PTempAll  = xyz_PTemp + xyz_PTempBZ
    xyz_ExnerAll  = xyz_Exner + xyz_ExnerBZ
    xyzf_QMixAll  = xyzf_QMix + xyzf_QMixBZ

    ! Perturbation component of Exner function at the surface is assumed 
    ! to be same as that at the lowest layer. (YOT, 2011/09/03)
    ! 
    ExnerBZSfc    = (PressSfc / PressBasis) ** (GasRDry / CpDry)
    xy_PressSfc   = (PressBasis * xyz_ExnerAll(:,:,kz))**(CpDry / GasRDry)

    ! xyz 格子点の速度の計算
    ! Calculate velocities at xyz grid points
    !
    xyz_VelX = xyz_avr_pyz(pyz_VelX)
    xyz_VelY = xyz_avr_xqz(xqz_VelY)

    ! 水平風速の絶対値の計算
    ! Calculate of absoluto horizontal velocities
    ! 
    xyz_AbsVel = SQRT( xyz_VelX**2 + xyz_VelY**2 + Vel0**2 )
    pyz_AbsVel = pyz_avr_xyz(xyz_AbsVel)
    xqz_AbsVel = xqz_avr_xyz(xyz_AbsVel)

    ! バルク $ R_i $ 数算出
    ! Calculate bulk $ R_i $
    !
    xy_SurfBulkRiNum = Grav / ( xyz_PTempAll(:,:,1) ) * ( xyz_PTempAll(:,:,1) -TempSfc / (ExnerBZSfc + xyz_Exner(:,:,1))) / max( xyz_AbsVel(:,:,1), VelMinForRi )**2 * z_dz(1) * 0.5d0

    ! バルク係数の計算
    ! Bulk coefficients are calculated.
    ! 
    call BulkCoef( xy_SurfBulkRiNum, xy_SurfRoughLength, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQmixBulkCoef )
    ! フラックスの計算
    ! Surface fluxes are calculated. 
    !
    xy_PTempFlux = - xy_SurfTempBulkCoef * xyz_AbsVel(:,:,kz) * ( xyz_PTempAll(:,:,kz) - TempSfc / ( ExnerBZSfc + xyz_Exner(:,:,kz) ) )

    xyf_QMixFlux = 0.0d0
    do s = 1, CondNum
      xyf_QMixFlux(:,:,IdxCG(s)) = - xy_SurfQmixBulkCoef * xyz_AbsVel(:,:,kz) * ( xyzf_QMixAll(:,:,kz,s) - SvapPress( SpcWetID(IdxCC(s)), TempSfc ) / ( xy_PressSfc ) * (MolWtWet(IdxCG(s)) / MolWtDry) )
    end do
    !
    xy_ExnerFlux = xy_DExnerDt_xy_xyf(xy_PTempFlux, xyf_QMixFlux, kz)
    !
    py_VelXFlux = - xy_SurfVelBulkCoef * pyz_AbsVel(:,:,kz) * pyz_VelX(:,:,kz)
    !
    xq_VelYFlux = - xy_SurfVelBulkCoef * xqz_AbsVel(:,:,kz) * xqz_VelY(:,:,kz)


    ! フラックスの下限値を設定
    ! Set lower limit of surface fluxes
    ! 
    xy_PTempFlux = max( PTempFluxMin, xy_PTempFlux )
    xy_ExnerFlux = max( ExnerFluxMin, xy_ExnerFlux )
    xyf_QMixFlux = max( QMixFluxMin, xyf_QMixFlux )


    ! 地表フラックスによる時間変化を計算
    ! Tendencies by surface fluxes (convergences of fluxes) are calculated.
    !
    xy_DPTempDtBulk = - ( 0.0d0 - xy_PTempFlux ) / z_dz(kz)
    !
    xy_DExnerDtBulk = - ( 0.0d0 - xy_ExnerFlux ) / z_dz(kz)
    !
    xyf_DQMixDtBulk = - ( 0.0d0 - xyf_QMixFlux ) / z_dz(kz)
    !
    py_DVelXDtBulk = - ( 0.0d0 - py_VelXFlux ) / z_dz(kz)
    !
    xq_DVelYDtBulk = - ( 0.0d0 - xq_VelYFlux ) / z_dz(kz)

    
    ! 仮引数配列へ格納
    ! Add tendency by surface flux convergence
    !
    xyz_DPTempDt(:,:,kz) = xyz_DPTempDt(:,:,kz) + xy_DPTempDtBulk
    xyz_DExnerDt(:,:,kz) = xyz_DExnerDt(:,:,kz) + xy_DExnerDtBulk
    do s = 1, ncmax
      xyzf_DQMixDt(:,:,kz,s) = xyzf_DQMixDt(:,:,kz,s) + xyf_DQMixDtBulk(:,:,s)
    end do
    pyz_DVelXDt (:,:,kz) = pyz_DVelXDt (:,:,kz) + py_DVelXDtBulk
    xqz_DVelYDt (:,:,kz) = xqz_DVelYDt (:,:,kz) + xq_DVelYDtBulk

    ! 出力
    ! Output
    !
    xyz_DPTempDtBulk = 0.0d0
    xyz_DExnerDtBulk = 0.0d0
    xyzf_DQMixDtBulk = 0.0d0
    pyz_DVelXDtBulk  = 0.0d0
    xqz_DVelYDtBulk  = 0.0d0

    xyz_DPTempDtBulk(:,:,kz) = xy_DPTempDtBulk
    xyz_DExnerDtBulk(:,:,kz) = xy_DExnerDtBulk
    do s = 1, ncmax
      xyzf_DQMixDtBulk(:,:,kz,s) = xyf_DQMixDtBulk(:,:,s)
    end do
    pyz_DVelXDtBulk (:,:,kz) = py_DVelXDtBulk
    xqz_DVelYDtBulk (:,:,kz) = xq_DVelYDtBulk
    !
    call HistoryAutoPut(TimeN, 'PTempSfc', xyz_DPTempDtBulk(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerSfc', xyz_DExnerDtBulk(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSfc',  pyz_DVelXDtBulk (1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSfc',  xqz_DVelYDtBulk (1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Sfc', xyzf_DQMixDtBulk(1:nx,1:ny,1:nz,s))
    end do


    call HistoryAutoPut(TimeN, 'PTempSfcFlux', xy_PTempFlux(1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'ExnerSfcFlux', xy_ExnerFlux(1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'VelXSfcFlux',  py_VelXFlux (1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'VelYSfcFlux',  xq_VelYFlux (1:nx,1:ny))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcFlux', xyf_QMixFlux(1:nx,1:ny,s))
    end do

    call HistoryAutoPut(TimeN, 'SfcHeatFlux', CpDry * xyz_DensBZ(1:nx,1:ny,1) * xy_PTempFlux(1:nx,1:ny) * ( ExnerBZSfc + xyz_Exner(1:nx,1:ny,1) ) )
    call HistoryAutoPut(TimeN, 'SfcXMomFlux', xyz_DensBZ(1:nx,1:ny,1) * py_VelXFlux (1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'SfcYMomFlux', xyz_DensBZ(1:nx,1:ny,1) * xq_VelYFlux (1:nx,1:ny))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcMassFlux', xyz_DensBZ(1:nx,1:ny,1) * xyf_QMixFlux(1:nx,1:ny,s))
    end do

  end subroutine Surfaceflux_Bulk_forcing

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

    ! 暗黙の型宣言禁止
    ! Implicit none
    !
    implicit none

    ! 作業変数
    ! Work variables
    !
    integer    :: l, unit   ! 出力装置番号
                         ! Device number 

    !---------------------------------------------------------------
    ! NAMELIST から情報を取得
    !
    NAMELIST /surfaceflux_bulk_nml/ FlagConstBulkCoef, FlagUseOfBulkCoefInNeutralCond, ConstBulkCoef, VelMinForRi, SfcRoughLength, Vel0, VelBulkCoefMin, TempBulkCoefMin, QmixBulkCoefMin, VelBulkCoefMax, TempBulkCoefMax, QmixBulkCoefMax             


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

    if (myrank == 0) then 
      call MessageNotify( "M", module_name, "SfcRoughLength = %f", d=(/SfcRoughLength/))
      call MessageNotify( "M", module_name, "VelMinForRi = %f", d=(/VelMinForRi/) )
      call MessageNotify( "M", module_name, "Vel0 = %f", d=(/Vel0/))

      call MessageNotify( 'M', module_name, "FlagConstBulkCoef              = %b", l = (/ FlagConstBulkCoef /) )
      call MessageNotify( 'M', module_name, "FlagUseOfBulkCoefInNeutralCond = %b", l = (/ FlagUseOfBulkCoefInNeutralCond /) )
      call MessageNotify( 'M', module_name, "ConstBulkCoef   = %f", d = (/ ConstBulkCoef   /) )
      call MessageNotify( 'M', module_name, "VelBulkCoefMin  = %f", d = (/ VelBulkCoefMin  /) )
      call MessageNotify( 'M', module_name, "TempBulkCoefMin = %f", d = (/ TempBulkCoefMin /) )
      call MessageNotify( 'M', module_name, "QmixBulkCoefMin = %f", d = (/ QmixBulkCoefMin /) )
      call MessageNotify( "M", module_name, "VelBulkCoefMax = %f", d=(/VelBulkCoefMax/))
      call MessageNotify( "M", module_name, "TempBulkCoefMax = %f", d=(/TempBulkCoefMax/))
      call MessageNotify( "M", module_name, "QmixBulkCoefMax = %f", d=(/QmixBulkCoefMax/))
    end if

    call HistoryAutoAddVariable( varname='PTempSfcFlux', dims=(/'x','y','t'/), longname='surface potential temperature flux (heat flux divided by density and specific heat)', units='K.m.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerSfcFlux', dims=(/'x','y','t'/), longname='surface exner function flux (heat flux divided by density and specific heat)', units='s-1', xtype='float')

    call HistoryAutoAddVariable( varname='VelXSfcFlux', dims=(/'x','y','t'/), longname='surface flux of x-component of velocity (momentum flux divided by density)', units='m2.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelYSfcFlux', dims=(/'x','y','t'/), longname='surface flux of y-component of velocity (momentum flux divided by density)', units='m2.s-2', xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_SfcFlux', dims=(/'x','y','t'/), longname='surface flux of ' //trim(SpcWetSymbol(l))//' mixing ratio (mass flux divided by density)', units='m.s-1', xtype='float')
    end do


    call HistoryAutoAddVariable( varname='PTempSfc', dims=(/'x','y','z','t'/), longname='potential temperature tendency by surface flux', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerSfc', dims=(/'x','y','z','t'/), longname='exner function tendency by surface flux', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='VelXSfc', dims=(/'x','y','z','t'/), longname='x-component velocity tendency by surface flux', units='m.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelYSfc', dims=(/'x','y','z','t'/), longname='y-component velocity tendency by surface flux', units='m.s-2', xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_Sfc', dims=(/'x','y','z','t'/), longname=trim(SpcWetSymbol(l))//' mixing ratio tendency by surface flux', units='s-1', xtype='float')
    end do


    call HistoryAutoAddVariable( varname='SfcHeatFlux', dims=(/'x','y','t'/), longname='surface heat flux', units='W.m-2', xtype='float')

    call HistoryAutoAddVariable( varname='SfcXMomFlux', dims=(/'x','y','t'/), longname='surface x-component momentum flux', units='kg.m-2.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='SfcYMomFlux', dims=(/'x','y','t'/), longname='surface y-component momentum flux', units='kg.m-2.s-1', xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_SfcMassFlux', dims=(/'x','y','t'/), longname=trim(SpcWetSymbol(l))//' surface mass flux', units='kg.m-2.s-1', xtype='float')
    end do

  end subroutine Surfaceflux_Bulk_init