!= 地表面フラックス (バルク法)
!
!= Surface flux (Bulk method)
!
! Authors::   Yasuhiro MORIKAWA, Yukiko YAMADA, Yoshiyuki O. Takahashi
! Version::   $Id: surface_flux_bulk.f90,v 1.27 2015/02/06 11:25:14 yot Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2009. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module surface_flux_bulk
  !
  != 地表面フラックス
  !
  != Surface flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 地表面フラックスを計算. 
  !
  ! Surface fluxes are calculated.
  !
  !== References
  !
  ! Louis, J-F., M. Tiedtke, and J-F. Geleyn, 
  ! A short history of the PBL parameterization at ECMWF, 
  ! Workshop on Planetary Boundary Layer Parameterization, 59-80, ECMWF, Reading, U.K., 
  ! 1982.
  !
  ! Beljaars, A. C. M., and A. A. M. Holtslag, 
  ! Flux parameterization over land surfaces for atmospheric models,
  ! J. Appl. Meteor., 30, 327-341, 1991.
  !
  ! Beljaars, A. C. M., 
  ! The parameterization of surface fluxes in large-scale models 
  ! under free convection, 
  ! Q. J. R. Meteorol. Soc., 121, 255-270, 1994.
  !
  !== Procedures List
  !
  ! SurfaceFlux       :: 地表面フラックスの計算
  ! SurfaceFluxOutput :: 地表面フラックスの出力
  ! ------------      :: ------------
  ! SurfaceFlux       :: Calculate surface fluxes
  ! SurfaceFluxOutput :: Output surface fluxes
  !
  !== NAMELIST
  !
  ! NAMELIST#surface_flux_bulk_nml
  !

  ! モジュール引用 ; USE statements
  !

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数. 
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数. 
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数. 
                             ! Number of vertical level

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  use composition, only: ncmax, IndexH2OVap

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING     ! 文字列.       Strings. 

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: SurfaceFlux
  public :: SurfaceFluxOutput
  public :: SurfaceFluxInit

  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  integer, parameter :: IDBulkCoefMethodL82     = 1
  integer, parameter :: IDBulkCoefMethodBH91B94 = 2

  logical            :: FlagIncludeB94W

  logical, save :: surface_flux_bulk_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  real(DP), save:: VelMinForRi
                            ! $ R_i $ 数用風最小値. 
                            ! Minimum value of velocity for $ R_i $ number
  real(DP), save:: VelMinForVel
                            ! 運動量用風最小値. 
                            ! Minimum value of velocity for momentum
  real(DP), save:: VelMinForTemp
                            ! 熱用風最小値. 
                            ! Minimum value of velocity for thermal
  real(DP), save:: VelMinForQVap
                            ! 水蒸気用風最小値. 
                            ! Minimum value of velocity for vapor
  real(DP), save:: VelMaxForVel
                            ! 運動量用風最大値. 
                            ! Maximum value of velocity for momentum
  real(DP), save:: VelMaxForTemp
                            ! 熱用風最大値. 
                            ! Maximum value of velocity for thermal
  real(DP), save:: VelMaxForQVap
                            ! 水蒸気用風最大値. 
                            ! Maximum value of velocity for vapor


  ! バルク係数
  ! Bluk coefficients
  !
  logical, save:: FlagConstBulkCoef
                            ! Flag for using constant bulk coefficient
  logical, save:: FlagUseOfBulkCoefInNeutralCond
                            ! Flag for using bulk coefficient in neutral condition
  real(DP), save:: ConstBulkCoef
                            ! バルク係数一定値. 
                            ! Steady value of bulk coefficient
  real(DP), save:: VelBulkCoefMin
                            ! $ u $ バルク係数最小値. 
                            ! Minimum value of $ u $ bulk coefficient
  real(DP), save:: TempBulkCoefMin
                            ! $ T $ バルク係数最小値. 
                            ! Minimum value of $ T $ bulk coefficient
  real(DP), save:: QVapBulkCoefMin
                            ! $ q $ バルク係数最小値. 
                            ! Minimum value of $ q $ bulk coefficient
  real(DP), save:: VelBulkCoefMax
                            ! $ u $ バルク係数最大値. 
                            ! Maximum value of $ u $ bulk coefficient
  real(DP), save:: TempBulkCoefMax
                            ! $ T $ バルク係数最大値. 
                            ! Maximum value of $ T $ bulk coefficient
  real(DP), save:: QVapBulkCoefMax
                            ! $ q $ バルク係数最大値. 
                            ! Maximum value of $ q $ bulk coefficient
  logical , save:: FlagFixFricTimeConstAtLB
  real(DP), save:: FricTimeConstAtLB
                            ! 下部境界摩擦の時定数 (s).
                            ! Time constant of surface friction (s).
  real(DP), save:: LowLatFricAtLB
                            ! 下部境界摩擦が働く最低緯度 (degree).
                            ! Lowest latitude where the friction is applied (degree)
  logical , save:: FlagFixHeatFluxAtLB
  real(DP), save:: HeatFluxAtLB
                            ! 下部境界での熱フラックス (W m-2).
                            ! Heat flux at the lower boundary (W m-2).
  logical , save:: FlagFixMassFluxAtLB
  real(DP), save:: MassFluxAtLB
                            ! 下部境界での質量フラックス (W m-2).
                            ! 実際にはゼロに固定するために使う程度にしか使えないだろう.
                            ! Mass flux at the lower boundary (kg m-2 s-1).

  character(*), parameter:: module_name = 'surface_flux_bulk'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: surface_flux_bulk.f90,v 1.27 2015/02/06 11:25:14 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version


contains

  !--------------------------------------------------------------------------------------

  subroutine SurfaceFlux( &
    & BulkCoefMethod,                                               & ! (in)
    & xyz_U, xyz_V,                                                 & ! (in)
    & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xy_SurfVirTemp, & ! (in)
    & xyzf_QMix,                                                    & ! (in)
    & xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner,   & ! (in)
    & xy_SurfTemp, xy_SurfHumidCoef,                                & ! (in)
    & xy_SurfRoughLengthMom, xy_SurfRoughLengthHeat,                & ! (in)
    & xy_SnowFrac,                                                  & ! (in)
    & xy_MomFluxX, xy_MomFluxY, xy_HeatFlux, xyf_QMixFlux,          & ! (out)
    & xy_SurfVelTransCoef, xy_SurfTempTransCoef,                    & ! (out)
    & xy_SurfQVapTransCoef,                                         & ! (out)
    & xy_SurfMOLength                                               & ! (out)
    & )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! Calculate radiation flux from temperature, specific humidity, and 
    ! air pressure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI
                              ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & GasRDry, &
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air
      & CpDry
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: &
      & xy_CalcQVapSatOnLiq, &
      & xy_CalcQVapSatOnSol

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & y_Lat                 ! $ \varphi $ [rad.] . 緯度. Latitude

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

    ! デバッグ用ユーティリティ
    ! Utilities for debug
    !
    use dc_trace, only: DbgMessage, BeginSub, EndSub

    ! 宣言文 ; Declaration statements
    !
    implicit none

    character(*), intent(in):: BulkCoefMethod
                              !
                              ! Method for calculating bulk coefficient
    real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u $ . 東西風速. Eastward wind
    real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v $ . 南北風速. Northward wind

    real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ . 温度 (整数レベル). 
                              ! Temperature (full level)
    real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ T $ . 温度 (半整数レベル). 
                              ! Temperature (half level)
    real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_v $ . 仮温度 (整数レベル). 
                              ! Virtual temperature (full level)
    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ T_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: xy_SurfVirTemp (0:imax-1, 1:jmax)
                              ! $ T_v $ . 仮温度 (惑星表面). 
                              ! Virtual temperature (surface)
    real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q $ .     比湿. Specific humidity
    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ p_s $ . 地表面気圧 (半整数レベル). 
                              ! Surface pressure (half level)
    real(DP), intent(in):: xy_SurfHeight(0:imax-1,1:jmax)
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP), intent(in):: xy_SurfRoughLengthMom (0:imax-1, 1:jmax)
                              ! 地表粗度長. 
                              ! Surface rough length for momentum
    real(DP), intent(in):: xy_SurfRoughLengthHeat(0:imax-1, 1:jmax)
                              ! 地表粗度長. 
                              ! Surface rough length for heat
    real(DP), intent(in):: xy_SnowFrac(0:imax-1, 1:jmax)
                              ! 
                              ! Snow fraction
    real(DP), intent(out):: xy_MomFluxX (0:imax-1, 1:jmax)
                              ! 惑星表面東西方向運動量フラックス. 
                              ! Eastward momentum flux at surface
    real(DP), intent(out):: xy_MomFluxY (0:imax-1, 1:jmax)
                              ! 惑星表面南北方向運動量フラックス. 
                              ! Northward momentum flux at surface
    real(DP), intent(out):: xy_HeatFlux (0:imax-1, 1:jmax)
                              ! 惑星表面熱フラックス. 
                              ! Heat flux at surface
    real(DP), intent(out):: xyf_QMixFlux(0:imax-1, 1:jmax, 1:ncmax)
                              ! 惑星表面比湿フラックス. 
                              ! Specific humidity flux at surface
    real(DP), intent(out):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：水蒸気
                              ! Transfer coefficient: water vapor
    real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)

    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_SurfTempBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：温度. 
                              ! Bulk coefficient: temperature
    real(DP):: xy_SurfQVapBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：比湿. 
                              ! Bulk coefficient: specific humidity
    real(DP):: xy_SurfVelBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：運動量. 
                              ! Bulk coefficient: temperature
    real(DP):: xy_SurfVelAbs (0:imax-1, 1:jmax)
                              ! 風速絶対値. 
                              ! Absolute velocity
    real(DP):: xy_SurfQVapSatOnLiq (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface

    real(DP):: xy_MomFluxXSurf (0:imax-1, 1:jmax)
                              ! 地表面の東西方向運動量フラックス. 
                              ! Eastward momentum flux on surface
    real(DP):: xy_MomFluxYSurf (0:imax-1, 1:jmax)
                              ! 地表面の南北方向運動量フラックス. 
                              ! Northward momentum flux on surface
    real(DP):: xy_HeatFluxSurf (0:imax-1, 1:jmax)
                              ! 地表面の熱フラックス. 
                              ! Heat flux on surface
    real(DP):: xyf_QMixFluxSurf(0:imax-1, 1:jmax, 1:ncmax)
                              ! 地表面の質量フラックス. 
                              ! Mass flux of constituents on surface

    real(DP):: xy_BetaW   (0:imax-1, 1:jmax)
                              ! 
                              ! "vertical velocity" (B94)

    integer            :: IDBulkCoefMethod

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. surface_flux_bulk_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! Check method for calculating bulk coefficient
    !
    if ( BulkCoefMethod == 'L82' ) then
      IDBulkCoefMethod = IDBulkCoefMethodL82
    else if ( BulkCoefMethod == 'BH91B94' ) then
      IDBulkCoefMethod = IDBulkCoefMethodBH91B94
    else
      call MessageNotify( 'E', module_name, 'BulkCoefMethod of %c is inappropriate.', c1 = trim( BulkCoefMethod ) )
    end if


    ! バルク係数算出
    ! Calculate bulk coefficients
    !
    call BulkCoef( &
      & IDBulkCoefMethod, &       ! (in)
      & xy_SurfRoughLengthMom ,  & ! (in)
      & xy_SurfRoughLengthHeat,  & ! (in)
      & xy_SurfHeight,       & ! (in)
      & xyz_Height,          & ! (in)
      & xyz_U(:,:,1), xyz_V(:,:,1), & ! (in)
!!$      & xy_SurfTemp, xyz_Temp(:,:,1), xyr_Exner(:,:,0), xyz_Exner(:,:,1), & ! (in)
      & xy_SurfVirTemp, xyz_VirTemp(:,:,1), xyr_Exner(:,:,0), xyz_Exner(:,:,1), & ! (in)
      & xy_SurfVelBulkCoef,  & ! (out)
      & xy_SurfTempBulkCoef, & ! (out)
      & xy_SurfQVapBulkCoef, & ! (out)
      & xy_BetaW,            & ! (out)
      & xy_SurfMOLength      & ! (out)
      & )


    ! 
    ! Calculation of wind speed
    !
    xy_SurfVelAbs = sqrt ( xyz_U(:,:,1)**2 + xyz_V(:,:,1)**2 + xy_BetaW**2 )
!!$    xy_SurfVelAbs = sqrt ( xyz_U(:,:,1)**2 + xyz_V(:,:,1)**2 )


    ! 輸送係数の計算
    ! Calculate transfer coefficient
    !
    if ( .not. FlagFixFricTimeConstAtLB ) then
      do i = 0, imax-1
        do j = 1, jmax
!!$          xy_SurfVelTransCoef(i,j) = &
!!$            &   xy_SurfVelBulkCoef(i,j) &
!!$            &   * xyr_Press(i,j,0) / ( GasRDry * xyr_Temp(i,j,0) ) &
!!$            &   * min( max( xy_SurfVelAbs(i,j), VelMinForVel ), VelMaxForVel )
          xy_SurfVelTransCoef(i,j) =                                  &
            &   xy_SurfVelBulkCoef(i,j)                               &
            &   * xyr_Press(i,j,0) / ( GasRDry * xyr_VirTemp(i,j,0) ) &
            &   * min( max( xy_SurfVelAbs(i,j), VelMinForVel ), VelMaxForVel )
        end do
      end do
    else
      do j = 1, jmax
        if ( abs( y_Lat(j) ) >= LowLatFricAtLB * PI / 180.0_DP ) then
          xy_SurfVelTransCoef(:,j) = 1.0_DP / FricTimeConstAtLB
        else
          xy_SurfVelTransCoef(:,j) = 0.0_DP
        end if
      end do
    end if

    if ( .not. FlagFixHeatFluxAtLB ) then
      do i = 0, imax-1
        do j = 1, jmax
!!$          xy_SurfTempTransCoef(i,j) = &
!!$            &   xy_SurfTempBulkCoef(i,j) &
!!$            &   * xyr_Press(i,j,0) / ( GasRDry * xyr_Temp(i,j,0) ) &
!!$            &   * min( max( xy_SurfVelAbs(i,j), VelMinForTemp ), VelMaxForTemp )
          xy_SurfTempTransCoef(i,j) =                                 &
            &   xy_SurfTempBulkCoef(i,j)                              &
            &   * xyr_Press(i,j,0) / ( GasRDry * xyr_VirTemp(i,j,0) ) &
            &   * min( max( xy_SurfVelAbs(i,j), VelMinForTemp ), VelMaxForTemp )
        end do
      end do
    else
      ! Set meaningless value.
      xy_SurfTempTransCoef = 0.0_DP
    end if

    if ( .not. FlagFixMassFluxAtLB ) then
      do i = 0, imax-1
        do j = 1, jmax
!!$          xy_SurfQVapTransCoef(i,j) =                                            &
!!$            &   xy_SurfQVapBulkCoef(i,j)                                         &
!!$            &   * xyr_Press(i,j,0) / ( GasRDry * xyr_Temp(i,j,0) )               &
!!$            &   * min( max( xy_SurfVelAbs(i,j), VelMinForQVap ), VelMaxForQVap )
          xy_SurfQVapTransCoef(i,j) =                                            &
            &   xy_SurfQVapBulkCoef(i,j)                                         &
            &   * xyr_Press(i,j,0) / ( GasRDry * xyr_VirTemp(i,j,0) )            &
            &   * min( max( xy_SurfVelAbs(i,j), VelMinForQVap ), VelMaxForQVap )
        end do
      end do
    else
      ! Set meaningless value.
      xy_SurfQVapTransCoef = 0.0_DP
    end if

    ! 飽和比湿の計算
    ! Calculate saturated specific humidity
    !
!!$    xy_SurfQVapSat = xy_CalcQVapSat( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSatOnLiq = xy_CalcQVapSatOnLiq( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSatOnSol = xy_CalcQVapSatOnSol( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSat = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq &
      & + xy_SnowFrac              * xy_SurfQVapSatOnSol


    ! 地表面フラックスの計算
    ! Calculate fluxes on flux
    !
    !   Momentum
    !
    xy_MomFluxXSurf = - xy_SurfVelTransCoef * xyz_U(:,:,1)
    xy_MomFluxYSurf = - xy_SurfVelTransCoef * xyz_V(:,:,1)

    !   Heat
    !
    if ( .not. FlagFixHeatFluxAtLB ) then
      xy_HeatFluxSurf = - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef &
        &                   * (   xyz_Temp(:,:,1) / xyz_Exner(:,:,1) &
        &                       - xy_SurfTemp     / xyr_Exner(:,:,0) )
    else
      xy_HeatFluxSurf = HeatFluxAtLB
    end if

    !   Mass
    !
    if ( .not. FlagFixMassFluxAtLB ) then
      xyf_QMixFluxSurf(:,:,IndexH2OVap) =                      &
        & - xy_SurfHumidCoef * xy_SurfQVapTransCoef(:,:)       &
        & * ( xyzf_QMix(:,:,1,IndexH2OVap) - xy_SurfQVapSat )
    else
      xyf_QMixFluxSurf(:,:,IndexH2OVap) = MassFluxAtLB
    end if
    !
    xyf_QMixFluxSurf(:,:,1:IndexH2OVap-1)     = 0.0_DP
    xyf_QMixFluxSurf(:,:,IndexH2OVap+1:ncmax) = 0.0_DP

    ! Surface flux of constituents except for water vapor is zero.
!!$    write( 6, * ) "MEMO: Surface flux of constituents except for water vapor is zero. (YOT, 2009/08/14)"


    ! フラックスの計算
    ! Calculate fluxes
    !
    xy_MomFluxX = xy_MomFluxXSurf
    xy_MomFluxY = xy_MomFluxYSurf
    xy_HeatFlux = xy_HeatFluxSurf
    do n = 1, ncmax
      xyf_QMixFlux(:,:,n) = xyf_QMixFluxSurf(:,:,n)
    end do

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'BulkCoefMom' , xy_SurfVelBulkCoef   )
    call HistoryAutoPut( TimeN, 'BulkCoefHeat', xy_SurfTempBulkCoef  )

    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine SurfaceFlux

  !--------------------------------------------------------------------------------------

  subroutine BulkCoef( &
    & IDBulkCoefMethod, &       ! (in)
    & xy_SurfRoughLengthMom ,  & ! (in)
    & xy_SurfRoughLengthHeat,  & ! (in)
    & xy_SurfHeight,          & ! (in)
    & xyz_Height,             & ! (in)
    & xy_U, xy_V,                                   & ! (in)
    & xy_SurfVirTemp, xy_VirTemp, xy_SurfExner, xy_Exner, & ! (in)
    & xy_SurfVelBulkCoef,     & ! (out)
    & xy_SurfTempBulkCoef,    & ! (out)
    & xy_SurfQVapBulkCoef,    & ! (out)
    & xy_BetaW,               & ! (out)
    & xy_SurfMOLength         & ! (out)
    & )
    !
    ! バルク係数を算出します.
    !
    ! Bulk coefficients are calculated.
    !

    ! モジュール引用 ; USE statements
    !

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & FKarm                 ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in):: IDBulkCoefMethod
                              !
                              !
    real(DP), intent(in):: xy_SurfRoughLengthMom (0:imax-1, 1:jmax)
                              ! 地表粗度長
                              ! Surface rough length for momentum
    real(DP), intent(in):: xy_SurfRoughLengthHeat(0:imax-1, 1:jmax)
                              ! 地表粗度長
                              ! Surface rough length for heat
    real(DP), intent(in):: xy_SurfHeight(0:imax-1,1:jmax)
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度. 
                              ! Height
    real(DP), intent(in):: xy_U (0:imax-1, 1:jmax)
                              !
                              ! Eastward wind velocity at lowest level
    real(DP), intent(in):: xy_V (0:imax-1, 1:jmax)
                              !
                              ! Northward wind velocity at lowest level
    real(DP), intent(in):: xy_SurfVirTemp (0:imax-1, 1:jmax)
                              !
                              ! Surface virtual temperature
    real(DP), intent(in):: xy_SurfExner(0:imax-1, 1:jmax)
                              !
                              ! Exner function at the surface
    real(DP), intent(in):: xy_VirTemp     (0:imax-1, 1:jmax)
                              !
                              ! Virtual temperature at lowest layer
    real(DP), intent(in):: xy_Exner    (0:imax-1, 1:jmax)
                              !
                              ! Exner function at lowest layer
    real(DP), intent(out):: xy_SurfVelBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：運動量. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfTempBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：温度. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfQVapBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：比湿. 
                              ! Bulk coefficient: specific humidity
    real(DP), intent(out):: xy_BetaW   (0:imax-1, 1:jmax)
                              ! 
                              ! "vertical velocity" (B94)
    real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)
                              !
                              ! Monin-Obukov length

    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_SurfBulkRiNum (0:imax-1, 1:jmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $ number
    real(DP):: xy_SurfVelAbs (0:imax-1, 1:jmax)
                              ! 風速絶対値. 
                              ! Absolute velocity
    real(DP) :: xy_SurfBulkCoefMomInNeutCond    (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfBulkCoefHeatInNeutCond   (0:imax-1, 1:jmax)

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude

    ! 実行文 ; Executable statement
    !

    if ( FlagConstBulkCoef ) then

      ! Use of constant bulk coefficient
      !

      xy_SurfVelBulkCoef  = ConstBulkCoef
      xy_SurfTempBulkCoef = ConstBulkCoef
      xy_SurfQVapBulkCoef = ConstBulkCoef

      xy_BetaW = 0.0_DP

      xy_SurfBulkRiNum = 0.0_DP

    else

      select case ( IDBulkCoefMethod )
      case ( IDBulkCoefMethodL82 )

        ! Parameterization by Louis et al. (1982)
        !

        ! 中立バルク係数の計算
        ! Calculate bulk coefficient in neutral condition
        !
        xy_SurfBulkCoefMomInNeutCond  =                     &
          & ( FKarm                                         &
          & / log (   ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthMom ) &
          &         / xy_SurfRoughLengthMom  ) )**2
        xy_SurfBulkCoefHeatInNeutCond  =                    &
          &   ( FKarm                                       &
          & / log (   ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthMom ) &
          &         / xy_SurfRoughLengthMom  ) )            &
          & * ( FKarm                                       &
          & / log (   ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthHeat ) &
          &         / xy_SurfRoughLengthHeat ) )

        if ( FlagUseOfBulkCoefInNeutralCond ) then

          ! 中立条件でのバルク係数の設定
          ! Set bulk coefficient in neutral condition
          !

          xy_SurfVelBulkCoef  = xy_SurfBulkCoefMomInNeutCond
          xy_SurfTempBulkCoef = xy_SurfBulkCoefHeatInNeutCond

          xy_SurfQVapBulkCoef = xy_SurfTempBulkCoef

          xy_SurfBulkRiNum = 0.0_DP

        else

          ! バルク $ R_i $ 数算出
          ! Calculate bulk $ R_i $
          !
          xy_SurfVelAbs = sqrt ( xy_U**2 + xy_V**2 )
          xy_SurfBulkRiNum =                                 &
            &   Grav / ( xy_SurfVirTemp / xy_SurfExner )     &
            &   * (   xy_VirTemp     / xy_Exner              &
            &       - xy_SurfVirTemp / xy_SurfExner )        &
            &   / max( xy_SurfVelAbs, VelMinForRi )**2       &
            &   * ( xyz_Height(:,:,1) - xy_SurfHeight )

          ! 非中立条件でのバルク係数の計算
          ! Calculate bulk coefficients in non-neutral condition
          !

          call BulkCoefL82( &
            & xy_SurfBulkRiNum, &       ! (in)
            & xy_SurfRoughLengthMom ,  & ! (in)
            & xy_SurfRoughLengthHeat,  & ! (in)
            & xy_SurfHeight,          & ! (in)
            & xyz_Height,             & ! (in)
            & xy_SurfBulkCoefMomInNeutCond,  & ! (in)
            & xy_SurfBulkCoefHeatInNeutCond, & ! (in)
            & xy_SurfVelBulkCoef,     & ! (out)
            & xy_SurfTempBulkCoef,    & ! (out)
            & xy_SurfQVapBulkCoef,    & ! (out)
            & xy_SurfMOLength         & ! (out)
            & )

        end if

        xy_BetaW = 0.0_DP

      case ( IDBulkCoefMethodBH91B94 )

        ! Parameterization by Beljaars and Holtslag (1991), Beljaars (1994)
        !

        call BulkCoefBH91B94( &
          & xy_SurfRoughLengthMom , & ! (in)
          & xy_SurfRoughLengthHeat, & ! (in)
          & xy_SurfHeight,          & ! (in)
          & xyz_Height,             & ! (in)
          & xy_U, xy_V,                                   & ! (in)
          & xy_SurfVirTemp, xy_VirTemp, xy_SurfExner, xy_Exner, & ! (in)
          & xy_SurfVelBulkCoef,     & ! (out)
          & xy_SurfTempBulkCoef,    & ! (out)
          & xy_SurfQVapBulkCoef,    & ! (out)
          & xy_BetaW,               & ! (out)
          & xy_SurfBulkRiNum,       & ! (out)
          & xy_SurfMOLength         & ! (out)
          & )

      end select

    end if


    ! 最大/最小 判定
    ! Measure maximum/minimum
    !
    do i = 0, imax-1
      do j = 1, jmax

        xy_SurfVelBulkCoef(i,j)  = &
          & max( min( xy_SurfVelBulkCoef(i,j), VelBulkCoefMax ), &
          &      VelBulkCoefMin )

        xy_SurfTempBulkCoef(i,j) = &
          & max( min( xy_SurfTempBulkCoef(i,j), TempBulkCoefMax ), &
          &      TempBulkCoefMin )

        xy_SurfQVapBulkCoef(i,j) = &
          & max( min( xy_SurfQVapBulkCoef(i,j), QVapBulkCoefMax ), &
          &      QVapBulkCoefMin )

      end do
    end do

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'SfcBulkRi', xy_SurfBulkRiNum )


  end subroutine BulkCoef

  !--------------------------------------------------------------------------------------

  subroutine BulkCoefL82( &
    & xy_SurfBulkRiNum, &       ! (in)
    & xy_SurfRoughLengthMom ,  & ! (in)
    & xy_SurfRoughLengthHeat,  & ! (in)
    & xy_SurfHeight,          & ! (in)
    & xyz_Height,             & ! (in)
    & xy_SurfBulkCoefMomInNeutCond, & ! (in)
    & xy_SurfBulkCoefHeatInNeutCond, & ! (in)
    & xy_SurfVelBulkCoef,     & ! (out)
    & xy_SurfTempBulkCoef,    & ! (out)
    & xy_SurfQVapBulkCoef,    & ! (out)
    & xy_SurfMOLength         & ! (out)
    & )
    !
    ! バルク係数を算出します.
    !
    ! Bulk coefficients are calculated.
    !

    ! モジュール引用 ; USE statements
    !
    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & FKarm                 ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in):: xy_SurfBulkRiNum (0:imax-1, 1:jmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $ number
    real(DP), intent(in):: xy_SurfRoughLengthMom (0:imax-1, 1:jmax)
                              ! 地表粗度長
                              ! Surface rough length for momentum
    real(DP), intent(in):: xy_SurfRoughLengthHeat(0:imax-1, 1:jmax)
                              ! 地表粗度長
                              ! Surface rough length for heat
    real(DP), intent(in):: xy_SurfHeight(0:imax-1,1:jmax)
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度. 
                              ! Height
    real(DP), intent(in):: xy_SurfBulkCoefMomInNeutCond (0:imax-1, 1:jmax)
                              !
                              !
    real(DP), intent(in):: xy_SurfBulkCoefHeatInNeutCond(0:imax-1, 1:jmax)
                              !
                              !
    real(DP), intent(out):: xy_SurfVelBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：運動量. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfTempBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：温度. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfQVapBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：比湿. 
                              ! Bulk coefficient: specific humidity
    real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)

    ! 作業変数
    ! Work variables
    !
    real(DP) :: SurfBulkRiNum
    real(DP) :: xy_MOLength(0:imax-1, 1:jmax)

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude

    ! 実行文 ; Executable statement
    !

    ! 非中立条件でのバルク係数の計算
    ! Calculate bulk coefficients in non-neutral condition
    !
    ! Parameterization by Louis et al. (1982)
    !
    do j = 1, jmax
      do i = 0, imax-1

        if ( xy_SurfBulkRiNum(i,j) > 0.0_DP ) then 

          xy_SurfVelBulkCoef(i,j) =                                       &
            &   xy_SurfBulkCoefMomInNeutCond(i,j)                         &
            &   / (   1.0_DP                                              &
            &       + 10.0_DP * xy_SurfBulkRiNum(i,j)                     &
            &         / sqrt( 1.0_DP + 5.0_DP * xy_SurfBulkRiNum(i,j) )   &
            &     )

          xy_SurfTempBulkCoef(i,j) =                                      &
            &   xy_SurfBulkCoefHeatInNeutCond(i,j)                        &
            &   / (   1.0_DP                                              &
            &       + 15.0_DP * xy_SurfBulkRiNum(i,j)                     &
            &         * sqrt( 1.0_DP + 5.0_DP * xy_SurfBulkRiNum(i,j) )   &
            &     )

          xy_SurfQVapBulkCoef(i,j) = xy_SurfTempBulkCoef(i,j)

        else

          xy_SurfVelBulkCoef(i,j) =                                              &
            &   xy_SurfBulkCoefMomInNeutCond(i,j)                                &
            &   * (   1.0_DP                                                     &
            &       - 10.0_DP * xy_SurfBulkRiNum(i,j)                            &
            &         / (   1.0_DP                                               &
            &             + 75.0_DP * xy_SurfBulkCoefMomInNeutCond(i,j)          &
            &               * sqrt( - ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) + xy_SurfRoughLengthMom(i,j) ) &
            &                         / xy_SurfRoughLengthMom(i,j)               &
            &                         * xy_SurfBulkRiNum(i,j)                    &
            &                     )                                              &
            &           )                                                        &
            &     )

          xy_SurfTempBulkCoef(i,j) =                                             &
            &   xy_SurfBulkCoefHeatInNeutCond(i,j)                               &
            &   * (   1.0_DP                                                     &
            &       - 15.0_DP * xy_SurfBulkRiNum(i,j)                            &
            &         / (   1.0_DP                                               &
            &             + 75.0_DP * xy_SurfBulkCoefHeatInNeutCond(i,j)         &
            &               * sqrt( - ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) + xy_SurfRoughLengthHeat(i,j) ) &
            &                         / xy_SurfRoughLengthHeat(i,j)              &
            &                         * xy_SurfBulkRiNum(i,j)                    &
            &                     )                                              &
            &           )                                                        &
            &     )

          xy_SurfQVapBulkCoef(i,j) = xy_SurfTempBulkCoef(i,j)

        end if

        ! Calculation of Monin-Obukhov length
        SurfBulkRiNum = xy_SurfBulkRiNum(i,j)
        if ( SurfBulkRiNum == 0.0_DP ) SurfBulkRiNum = 1.0e-10_DP
        xy_MOLength(i,j) =                               &
          &   ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) ) &
          & / ( FKarm * SurfBulkRiNum )                  &
          & * xy_SurfVelBulkCoef(i,j)**1.5_DP / xy_SurfTempBulkCoef(i,j)

      end do
    end do

    xy_SurfMOLength = xy_MOLength


  end subroutine BulkCoefL82

  !--------------------------------------------------------------------------------------

  subroutine BulkCoefBH91B94( &
    & xy_SurfRoughLengthMom ,                       & ! (in)
    & xy_SurfRoughLengthHeat,                       & ! (in)
    & xy_SurfHeight,                                & ! (in)
    & xyz_Height,                                   & ! (in)
    & xy_U, xy_V,                                   & ! (in)
    & xy_SurfVirTemp, xy_VirTemp, xy_SurfExner, xy_Exner, & ! (in)
    & xy_SurfVelBulkCoef,                           & ! (out)
    & xy_SurfTempBulkCoef,                          & ! (out)
    & xy_SurfQVapBulkCoef,                          & ! (out)
    & xy_BetaW,                                     & ! (out)
    & xy_SurfBulkRiNum,                             & ! (out)
    & xy_SurfMOLength                               & ! (out)
    & )
    !
    ! バルク係数を算出します.
    !
    ! Bulk coefficients are calculated.
    !

    ! モジュール引用 ; USE statements
    !

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

!!$    ! MPI 関連ルーチン
!!$    ! MPI related routines
!!$    !
!!$    use mpi_wrapper, only : MPIWrapperChkTrue

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & FKarm                 ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in):: xy_SurfRoughLengthMom (0:imax-1, 1:jmax)
                              ! 地表粗度長
                              ! Surface rough length for momentum
    real(DP), intent(in):: xy_SurfRoughLengthHeat(0:imax-1, 1:jmax)
                              ! 地表粗度長
                              ! Surface rough length for heat
    real(DP), intent(in):: xy_SurfHeight(0:imax-1,1:jmax)
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度. 
                              ! Height
    real(DP), intent(in):: xy_U (0:imax-1, 1:jmax)
                              !
                              ! Eastward wind velocity at lowest level
    real(DP), intent(in):: xy_V (0:imax-1, 1:jmax)
                              !
                              ! Northward wind velocity at lowest level
    real(DP), intent(in):: xy_SurfVirTemp (0:imax-1, 1:jmax)
                              !
                              ! Surface virtual temperature
    real(DP), intent(in):: xy_SurfExner(0:imax-1, 1:jmax)
                              !
                              ! Exner function at the surface
    real(DP), intent(in):: xy_VirTemp     (0:imax-1, 1:jmax)
                              !
                              ! Virtual temperature at lowest layer
    real(DP), intent(in):: xy_Exner    (0:imax-1, 1:jmax)
                              !
                              ! Exner function at lowest layer

    real(DP), intent(out):: xy_SurfVelBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：運動量. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfTempBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：温度. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfQVapBulkCoef (0:imax-1, 1:jmax)
                              ! バルク係数：比湿. 
                              ! Bulk coefficient: specific humidity
    real(DP), intent(out):: xy_BetaW   (0:imax-1, 1:jmax)
                              ! 
                              ! "vertical velocity" (B94)
    real(DP), intent(out):: xy_SurfBulkRiNum (0:imax-1, 1:jmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $ number
    real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)

    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_SurfVelAbs (0:imax-1, 1:jmax)
                              ! 風速絶対値. 
                              ! Absolute velocity
    real(DP) :: xy_MOLength    (0:imax-1, 1:jmax)
    real(DP) :: xy_MOLengthSave(0:imax-1, 1:jmax)
    real(DP) :: xy_ZetaM(0:imax-1, 1:jmax)
    real(DP) :: xy_ZetaH(0:imax-1, 1:jmax)
    real(DP) :: xy_PsiM1(0:imax-1, 1:jmax)
    real(DP) :: xy_PsiH1(0:imax-1, 1:jmax)
    real(DP) :: xy_PsiM0(0:imax-1, 1:jmax)
    real(DP) :: xy_PsiH0(0:imax-1, 1:jmax)
    real(DP) :: xy_FricVelByU1(0:imax-1, 1:jmax)
    real(DP) :: SurfBulkRiNum
    logical  :: FlagConverge
    logical  :: xy_FlagConverge(0:imax-1, 1:jmax)
    integer             :: iLoop
    integer , parameter :: nLoop = 100
    real(DP)            :: MOLengthErr
    real(DP), parameter :: MOLengthErrCriterion = 1.0d-5

    real(DP), parameter :: Gamma = 16.0_DP
    real(DP), parameter :: Beta  = 1.2_DP

    real(DP) :: SurfPotTemp
    real(DP) :: PotTemp
    real(DP) :: BLHeight
    real(DP) :: xy_BLHeight(0:imax-1, 1:jmax)

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude

    ! 実行文 ; Executable statement
    !

    ! Calculate bulk coefficients 
    ! Parameterization by Beljaars and Holtslag (1991)
    !

    ! initialization
    xy_MOLength = 1.0e10_DP
    xy_BetaW    = 0.0_DP

    xy_FlagConverge = .false.

    loop_iteration : do iLoop = 1, nLoop

      xy_MOLengthSave = xy_MOLength

      ! Calculation of Psi_{M,H}
      if ( FlagUseOfBulkCoefInNeutralCond ) then
        xy_PsiM0 = 0.0_DP
        xy_PsiH0 = 0.0_DP
        xy_PsiM1 = 0.0_DP
        xy_PsiH1 = 0.0_DP
      else
        xy_ZetaM = xy_SurfRoughLengthMom  / xy_MOLength
        xy_ZetaH = xy_SurfRoughLengthHeat / xy_MOLength
        call BH91CalcPsi(       &
          & Gamma,              & ! (in)
          & xy_ZetaM, xy_ZetaH, & ! (in)
          & xy_PsiM0, xy_PsiH0  & ! (out)
          & )
        xy_ZetaM = ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthMom  ) / xy_MOLength
        xy_ZetaH = ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthHeat ) / xy_MOLength
        call BH91CalcPsi(       &
          & Gamma,              & ! (in)
          & xy_ZetaM, xy_ZetaH, & ! (in)
          & xy_PsiM1, xy_PsiH1  & ! (out)
          & )
      end if

      do j = 1, jmax
        do i = 0, imax-1
          if ( .not. xy_FlagConverge(i,j) ) then
            ! u_* / U_1, Eq. (5)
            xy_FricVelByU1(i,j) =                             &
              &   FKarm                                       &
              & / (                                           &
              &       log (   ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) + xy_SurfRoughLengthMom(i,j)  ) &
              &             / xy_SurfRoughLengthMom(i,j)  )   &
              &     - xy_PsiM1(i,j) + xy_PsiM0(i,j)           &
              &   )

            xy_SurfVelBulkCoef(i,j) = xy_FricVelByU1(i,j)**2
            !
            xy_SurfTempBulkCoef(i,j) =                        &
              &   xy_FricVelByU1(i,j)                         &
              & * FKarm                                       &
              & / (                                           &
              &       log (   ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) + xy_SurfRoughLengthHeat(i,j) ) &
              &             / xy_SurfRoughLengthHeat(i,j) )   &
              &     - xy_PsiH1(i,j) + xy_PsiH0(i,j)           &
              &   )
            !
            xy_SurfQVapBulkCoef(i,j) = xy_SurfTempBulkCoef(i,j)
          end if
        end do
      end do


      ! Calculation of wind speed related to convection
      if ( FlagIncludeB94W ) then
        do j = 1, jmax
          do i = 0, imax-1
            if ( .not. xy_FlagConverge(i,j) ) then
              if ( xy_MOLength(i,j) < 0.0_DP ) then
                SurfPotTemp = xy_SurfVirTemp(i,j) / xy_SurfExner(i,j)
                PotTemp     = xy_VirTemp    (i,j) / xy_Exner    (i,j)
!!$              BLHeight =                                                      &
!!$                & - xy_MOLength(i,j) / ( FKarm**2 * Beta**3 )                 &
!!$                &     * (   log( - 38.5_DP * xy_MOLength(i,j)                 &
!!$                &                  / ( Gamma * xy_SurfRoughLengthMom(i,j) ) ) &
!!$                &         + xy_PsiM0(i,j) )**3
                ! BLHeight is assumed to be constant.
                BLHeight = 1000.0_DP
                xy_BetaW(i,j) =                                                 &
                  & sqrt(                                                       &
                  &       Grav / PotTemp * ( SurfPotTemp - PotTemp )            &
                  &         * FKarm**2 * Beta * BLHeight                        &
                  &     / (                                                     &
                  &           (   log( - 38.5_DP * xy_MOLength(i,j)             &
                  &                / ( Gamma * xy_SurfRoughLengthMom (i,j) ) )  &
                  &             + xy_PsiM0(i,j)                                )&
                  &         * (   log( -  4.0_DP * xy_MOLength(i,j)             &
                  &                / ( Gamma * xy_SurfRoughLengthHeat(i,j) ) )  &
                  &             + xy_PsiH0(i,j)                                )&
                  &       )                                                     &
                  &     )
                xy_BLHeight(i,j) = BLHeight
              else
                xy_BetaW(i,j) = 0.0_DP
                xy_BLHeight(i,j) = 0.0_DP
              end if
            end if
          end do
        end do
      else
        do j = 1, jmax
          do i = 0, imax-1
            if ( .not. xy_FlagConverge(i,j) ) then
              xy_BetaW(i,j) = 0.0_DP
            end if
          end do
        end do
      end if


      do j = 1, jmax
        do i = 0, imax-1
          if ( .not. xy_FlagConverge(i,j) ) then

            ! Calculation of bulk Richardson number
            xy_SurfVelAbs(i,j) = sqrt ( xy_U(i,j)**2 + xy_V(i,j)**2 + xy_BetaW(i,j)**2 )
            xy_SurfBulkRiNum(i,j) =                                  &
              &   Grav / ( xy_SurfVirTemp(i,j) / xy_SurfExner(i,j) ) &
              &   * (   xy_VirTemp(i,j)     / xy_Exner(i,j)          &
              &       - xy_SurfVirTemp(i,j) / xy_SurfExner(i,j) )    &
              &   / max( xy_SurfVelAbs(i,j), VelMinForRi )**2        &
              &   * ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) )

            ! Calculation of Monin-Obukhov length
            SurfBulkRiNum = xy_SurfBulkRiNum(i,j)
            if ( SurfBulkRiNum == 0.0_DP ) SurfBulkRiNum = 1.0e-10_DP
            xy_MOLength(i,j) =                               &
              &   ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) ) &
              & / ( FKarm * SurfBulkRiNum )                  &
              & * xy_SurfVelBulkCoef(i,j)**1.5_DP / xy_SurfTempBulkCoef(i,j)

          end if
        end do
      end do


      ! TO BE DELETED
!!$      ! Check of convergence
!!$      FlagConverge = .true.
!!$      loop_check : do j = 1, jmax
!!$        do i = 0, imax-1
!!$          MOLengthErr =                                      &
!!$            & abs( xy_MOLength(i,j) - xy_MOLengthSave(i,j) ) &
!!$            &    / max( abs( xy_MOLength(i,j) ), 1.0d-10 )
!!$          if ( MOLengthErr > MOLengthErrCriterion ) then
!!$            FlagConverge = .false.
!!$            exit loop_check
!!$          end if
!!$        end do
!!$      end do loop_check
!!$      a_FlagReCalcLocal = ( .not. FlagConverge )
!!$      call MPIWrapperChkTrue(   &
!!$        & 1, a_FlagReCalcLocal, & ! (in)
!!$        & a_FlagReCalcGlobal    & ! (out)
!!$        & )
!!$      if ( .not. a_FlagReCalcGlobal(1) ) exit loop_iteration

      ! Check of convergence
      do j = 1, jmax
        do i = 0, imax-1
          if ( .not. xy_FlagConverge(i,j) ) then
            MOLengthErr =                                      &
              & abs( xy_MOLength(i,j) - xy_MOLengthSave(i,j) ) &
              &    / max( abs( xy_MOLength(i,j) ), 1.0e-10_DP )
            if ( MOLengthErr <= MOLengthErrCriterion ) then
              xy_FlagConverge(i,j) = .true.
            end if
          end if
        end do
      end do
      FlagConverge = .true.
      loop_check : do j = 1, jmax
        do i = 0, imax-1
          if ( .not. xy_FlagConverge(i,j) ) then
            FlagConverge = .false.
            exit loop_check
          end if
        end do
      end do loop_check
      if ( FlagConverge ) exit loop_iteration


    end do loop_iteration
    if ( iLoop > nLoop ) then
      call MessageNotify( 'E', module_name, 'Monin-Obukhov length is not convergent.' )
    end if


    xy_SurfMOLength = xy_MOLength


    call HistoryAutoPut( TimeN, 'MOLength'   , xy_MOLength        )
    call HistoryAutoPut( TimeN, 'MOLengthInv', 1.0_DP/xy_MOLength )
    call HistoryAutoPut( TimeN, 'BetaW'      , xy_BetaW           )
    call HistoryAutoPut( TimeN, 'BLHeight'   , xy_BLHeight        )


  end subroutine BulkCoefBH91B94

  !--------------------------------------------------------------------------------------

  subroutine BH91CalcPsi( &
    & Gamma,              & ! (in)
    & xy_ZetaM, xy_ZetaH, & ! (in)
    & xy_PsiM, xy_PsiH    & ! (out)
    & )
    !
    ! 
    !
    ! Calculation of Psi_M and Psi_H
    !

    ! モジュール引用 ; USE statements
    !

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI
                              ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: Gamma

    real(DP), intent(in ) :: xy_ZetaM(0:imax-1, 1:jmax)
                              ! 
                              ! zeta = z / L for momentum
    real(DP), intent(in ) :: xy_ZetaH(0:imax-1, 1:jmax)
                              ! 
                              ! zeta = z / L for heat
    real(DP), intent(out) :: xy_PsiM(0:imax-1, 1:jmax)
                              ! 
                              ! PsiM
    real(DP), intent(out) :: xy_PsiH(0:imax-1, 1:jmax)
                              ! 
                              ! PsiH

    ! 作業変数
    ! Work variables
    !
    real(DP) :: ZetaM
    real(DP) :: ZetaH
    real(DP) :: ParamXM
    real(DP) :: ParamXH

    real(DP), parameter :: ConstA = 1.0_DP
    real(DP), parameter :: ConstB = 0.667_DP
    real(DP), parameter :: ConstC = 5.0_DP
    real(DP), parameter :: ConstD = 0.35_DP

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude

    ! 実行文 ; Executable statement
    !

    ! Parameterization by Beljaars and Holtslag (1991)
    !


    do j = 1, jmax
      do i = 0, imax-1

        ZetaM = xy_ZetaM(i,j)
        ZetaH = xy_ZetaH(i,j)

        if ( ZetaM < 0.0_DP ) then
          ! for unstable condition

          ParamXM = ( 1.0_DP - Gamma * ZetaM )**0.25_DP
          ParamXH = ( 1.0_DP - Gamma * ZetaH )**0.25_DP

          ! Eq. (25)
          xy_PsiM(i,j) =                                    &
            &   log(                                        &
            &          ( 1.0_DP + ParamXM    )**2           &
            &        * ( 1.0_DP + ParamXM**2 )    / 8.0_DP  &
            &      )                                        &
            & - 2.0_DP * atan( ParamXM ) + PI / 2.0_DP
          ! Eq. (26)
          xy_PsiH(i,j) =                                    &
            &   log(                                        &
            &          ( 1.0_DP + ParamXH**2 )**2 / 4.0_DP  &
            &      )
        else
          ! for stable condition

          ! Eq. (28)
          xy_PsiM(i,j) =                                                      &
            & - ConstA * ZetaM                                                &
            & - ConstB * ( ZetaM - ConstC / ConstD ) * exp( - ConstD * ZetaM )&
            & - ConstB * ConstC / ConstD
          ! Eq. (32)
          xy_PsiH(i,j) =                                                      &
            & - ( 1.0_DP + 2.0_DP / 3.0_DP * ConstA * ZetaH )**1.5            &
            & - ConstB * ( ZetaH - ConstC / ConstD ) * exp( - ConstD * ZetaH )&
            & - ConstB * ConstC / ConstD                                      &
            & + 1.0_DP
        end if

      end do
    end do


  end subroutine BH91CalcPsi

  !--------------------------------------------------------------------------------------

  subroutine SurfaceFluxOutput(                                &
    & xy_SnowFrac,                                             & ! (in)
    & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in)
    & xy_SurfH2OVapFluxA, xy_SurfLatentHeatFluxA,              & ! (in)
    & xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt,           & ! (in)
    & xy_SurfTemp, xy_DSurfTempDt,                             & ! (in)
    & xyr_Press, xyz_Exner, xyr_Exner, xy_SurfHumidCoef,       & ! (in)
    & xy_SurfVelTransCoef, xy_SurfTempTransCoef,               & ! (in)
    & xy_SurfQVapTransCoef                                     & ! (in)
    & )
    !
    ! フラックス (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux). 
    ! について, その他の引数を用いて補正し, 出力を行う. 
    !
    ! Fluxes (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux) are
    ! corrected by using other arguments, and the corrected values are output.
    !

    ! モジュール引用 ; USE statements
    !

    ! 物理定数設定
    ! Physical constant settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation

    ! 飽和比湿の算出
    ! Evaluation of saturation specific humidity
    !
    use saturate, only: &
      & xy_CalcQVapSatOnLiq,       &
      & xy_CalcQVapSatOnSol,       &
      & xy_CalcDQVapSatDTempOnLiq, &
      & xy_CalcDQVapSatDTempOnSol

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $ [s]
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xy_SnowFrac(0:imax-1, 1:jmax)
                              ! 
                              ! Snow fraction
    real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(in):: xy_SurfH2OVapFluxA    (0:imax-1, 1:jmax)
                              ! 惑星表面水蒸気フラックス.
                              ! Water vapor flux at the surface
    real(DP), intent(in):: xy_SurfLatentHeatFluxA(0:imax-1, 1:jmax)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface
    real(DP), intent(in):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{u}{t} $ . 東西風速時間変化率. 
                              ! Eastward wind tendency
    real(DP), intent(in):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{v}{t} $ . 南北風速時間変化率. 
                              ! Northward wind tendency
    real(DP), intent(in):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度時間変化率. 
                              ! Temperature tendency
    real(DP), intent(in):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 比湿時間変化率. 
                              ! Specific humidity tendency
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xy_DSurfTempDt (0:imax-1, 1:jmax)
                              ! 地表面温度時間変化率. 
                              ! Surface temperature tendency
    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)
    real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
                              ! 輸送係数：水蒸気
                              ! Transfer coefficient: water vapor

    ! 出力のための作業変数
    ! Work variables for output
    !
    real(DP):: xyr_MomFluxXCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP):: xyr_MomFluxYCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP):: xyr_HeatFluxCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP):: xyrf_QMixFluxCor(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP):: xyr_LatentHeatFluxCor(0:imax-1, 1:jmax, 0:kmax)
                              ! 表面潜熱フラックス.
                              ! Latent heat flux
    real(DP):: xy_SurfQVapSatOnLiq (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface
    real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface

    ! 作業変数
    ! Work variables
    !
    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. surface_flux_bulk_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! 飽和比湿の計算
    ! Calculate saturated specific humidity
    !
    xy_SurfQVapSatOnLiq  = &
      & xy_CalcQVapSatOnLiq( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSatOnSol  = &
      & xy_CalcQVapSatOnSol( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfQVapSat       = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq &
      & + xy_SnowFrac              * xy_SurfQVapSatOnSol
    xy_SurfDQVapSatDTempOnLiq = &
      & xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
    xy_SurfDQVapSatDTempOnSol = &
      & xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
    xy_SurfDQVapSatDTemp = &
      &   ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq &
      & + xy_SnowFrac              * xy_SurfDQVapSatDTempOnSol

    ! Output of fluxes at t
    !

    ! 風速, 温度, 比湿フラックス補正
    ! Correct fluxes of wind, temperature, specific humidity
    !
    do j = 1, jmax
      do i = 0, imax-1
        xyr_MomFluxXCor( i,j,0 ) = xyr_MomFluxX( i,j,0 ) &
          & - xy_SurfVelTransCoef( i,j ) * xyz_DUDt( i,j,1 ) * DelTime

        xyr_MomFluxYCor( i,j,0 ) = xyr_MomFluxY( i,j,0 ) &
          & - xy_SurfVelTransCoef( i,j ) * xyz_DVDt( i,j,1 ) * DelTime

        xyr_HeatFluxCor( i,j,0 ) = xyr_HeatFlux( i,j,0 )               &
          & - CpDry * xyr_Exner( i,j,0 ) * xy_SurfTempTransCoef( i,j ) &
          &     * ( xyz_DTempDt( i,j,1 ) / xyz_Exner( i,j,1 )          &
          &       - xy_DSurfTempDt( i,j ) / xyr_Exner( i,j,0 ) )       &
          &   * DelTime
      end do
    end do
    n = IndexH2OVap
    do j = 1, jmax
      do i = 0, imax-1
        xyrf_QMixFluxCor( i,j,0,n ) = xyrf_QMixFlux( i,j,0,n )                    &
          & - xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j )               &
          &   * ( xyzf_DQMixDt( i,j,1,n )                                         &
          &     - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * DelTime
      end do
    end do
    do n = 1, IndexH2OVap-1
      xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
    end do
    do n = IndexH2OVap+1, ncmax
      xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
    end do
    n = IndexH2OVap
    do j = 1, jmax
      do i = 0, imax-1
!!$        xyr_LatentHeatFluxCor( i,j,0 ) = xy_SurfLatentHeatFlux( i,j )             &
!!$          & - LatentHeat * xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j )  &
!!$          &   * ( xyzf_DQMixDt( i,j,1,n )                                         &
!!$          &     - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * DelTime
        xyr_LatentHeatFluxCor( i,j,0 ) = LatentHeat * xyrf_QMixFluxCor( i,j,0,n )
      end do
    end do


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'TauX'          , xyr_MomFluxXCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'TauY'          , xyr_MomFluxYCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'Sens'          , xyr_HeatFluxCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'SurfH2OVapFlux', xyrf_QMixFluxCor(:,:,0,IndexH2OVap) )
    call HistoryAutoPut( TimeN, 'Evap'          , xyr_LatentHeatFluxCor(:,:,0) )


    ! Output of fluxes at t - \Delta t
    !

    ! 風速, 温度, 比湿フラックス補正
    ! Correct fluxes of wind, temperature, specific humidity
    !
    do j = 1, jmax
      do i = 0, imax-1
        xyr_MomFluxXCor( i,j,0 ) = xyr_MomFluxX( i,j,0 )
        xyr_MomFluxYCor( i,j,0 ) = xyr_MomFluxY( i,j,0 )
        xyr_HeatFluxCor( i,j,0 ) = xyr_HeatFlux( i,j,0 )
      end do
    end do
    n = IndexH2OVap
    do j = 1, jmax
      do i = 0, imax-1
        xyrf_QMixFluxCor( i,j,0,n ) = xyrf_QMixFlux( i,j,0,n )
      end do
    end do
    do n = 1, IndexH2OVap-1
      xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
    end do
    do n = IndexH2OVap+1, ncmax
      xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
    end do
    n = IndexH2OVap
    do j = 1, jmax
      do i = 0, imax-1
!!$        xyr_LatentHeatFluxCor( i,j,0 ) = xy_SurfLatentHeatFlux( i,j )
        xyr_LatentHeatFluxCor( i,j,0 ) = LatentHeat * xyrf_QMixFluxCor( i,j,0,n )
      end do
    end do

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'TauXB'          , xyr_MomFluxXCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'TauYB'          , xyr_MomFluxYCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'SensB'          , xyr_HeatFluxCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'SurfH2OVapFluxB', xyrf_QMixFluxCor(:,:,0,IndexH2OVap) )
    call HistoryAutoPut( TimeN, 'EvapB'          , xyr_LatentHeatFluxCor(:,:,0) )


    ! Output of fluxes at t + \Delta t
    !

    ! 風速, 温度, 比湿フラックス補正
    ! Correct fluxes of wind, temperature, specific humidity
    !
    do j = 1, jmax
      do i = 0, imax-1
        xyr_MomFluxXCor( i,j,0 ) = xyr_MomFluxX( i,j,0 ) &
          & - xy_SurfVelTransCoef( i,j ) * xyz_DUDt( i,j,1 ) * 2.0_DP * DelTime

        xyr_MomFluxYCor( i,j,0 ) = xyr_MomFluxY( i,j,0 ) &
          & - xy_SurfVelTransCoef( i,j ) * xyz_DVDt( i,j,1 ) * 2.0_DP * DelTime

        xyr_HeatFluxCor( i,j,0 ) = xyr_HeatFlux( i,j,0 )               &
          & - CpDry * xyr_Exner( i,j,0 ) * xy_SurfTempTransCoef( i,j ) &
          &     * ( xyz_DTempDt( i,j,1 ) / xyz_Exner( i,j,1 )          &
          &       - xy_DSurfTempDt( i,j ) / xyr_Exner( i,j,0 ) )       &
          &   * 2.0_DP * DelTime
      end do
    end do
    n = IndexH2OVap
    do j = 1, jmax
      do i = 0, imax-1
        xyrf_QMixFluxCor( i,j,0,n ) = xyrf_QMixFlux( i,j,0,n )                    &
          & - xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j )               &
          &   * ( xyzf_DQMixDt( i,j,1,n )                                         &
          &     - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * 2.0_DP * DelTime
      end do
    end do
    do n = 1, IndexH2OVap-1
      xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
    end do
    do n = IndexH2OVap+1, ncmax
      xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
    end do
    n = IndexH2OVap
    do j = 1, jmax
      do i = 0, imax-1
!!$        xyr_LatentHeatFluxCor( i,j,0 ) = xy_SurfLatentHeatFlux( i,j )             &
!!$          & - LatentHeat * xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j )  &
!!$          &   * ( xyzf_DQMixDt( i,j,1,n )                                         &
!!$          &     - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * 2.0d0 * DelTime
        xyr_LatentHeatFluxCor( i,j,0 ) = LatentHeat * xyrf_QMixFluxCor( i,j,0,n )
      end do
    end do

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'TauXA'          , xyr_MomFluxXCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'TauYA'          , xyr_MomFluxYCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'SensA'          , xyr_HeatFluxCor (:,:,0) )
    call HistoryAutoPut( TimeN, 'SurfH2OVapFluxA', xyrf_QMixFluxCor(:,:,0,IndexH2OVap) )
    call HistoryAutoPut( TimeN, 'EvapA'          , xyr_LatentHeatFluxCor(:,:,0) )


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'SurfH2OVapFluxU', xy_SurfH2OVapFluxA     )
    call HistoryAutoPut( TimeN, 'EvapU'          , xy_SurfLatentHeatFluxA )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine SurfaceFluxOutput

  !--------------------------------------------------------------------------------------

  subroutine SurfaceFluxInit
    !
    ! surface_flux_bulk モジュールの初期化を行います. 
    ! NAMELIST#surface_flux_bulk_nml の読み込みはこの手続きで行われます. 
    !
    ! "surface_flux_bulk" module is initialized. 
    ! "NAMELIST#surface_flux_bulk_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: StoA

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & AxnameX, &
      & AxnameY, &
      & AxnameZ, &
      & AxnameR, &
      & AxnameT

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: SaturateInit


    ! 宣言文 ; Declaration statements
    !
    implicit none

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /surface_flux_bulk_nml/                                  &
      & FlagConstBulkCoef,                                            &
      & FlagUseOfBulkCoefInNeutralCond, ConstBulkCoef,                &
      & FlagIncludeB94W,                                              &
      !
      & VelMinForRi, VelMinForVel, VelMinForTemp, VelMinForQVap,      &
      & VelMaxForVel, VelMaxForTemp, VelMaxForQVap,                   &
      !
      & VelBulkCoefMin, TempBulkCoefMin, QVapBulkCoefMin,             &
      & VelBulkCoefMax, TempBulkCoefMax, QVapBulkCoefMax,             &
      !
      & FlagFixFricTimeConstAtLB, FricTimeConstAtLB, LowLatFricAtLB,  &
      & FlagFixHeatFluxAtLB,      HeatFluxAtLB,                       &
      & FlagFixMassFluxAtLB,      MassFluxAtLB
          !
          ! デフォルト値については初期化手続 "surface_flux_bulk#SurfaceFluxInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "surface_flux_bulk#SurfaceFluxInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( surface_flux_bulk_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    FlagConstBulkCoef              = .false.
    FlagUseOfBulkCoefInNeutralCond = .false.
    ConstBulkCoef                  =  0.0_DP
!!$    FlagIncludeB94W                = .false.
    FlagIncludeB94W                = .true.

    VelMinForRi   = 0.01_DP
    VelMinForVel  = 0.01_DP
    VelMinForTemp = 0.01_DP
    VelMinForQVap = 0.01_DP
    VelMaxForVel  = 1000.0_DP
    VelMaxForTemp = 1000.0_DP
    VelMaxForQVap = 1000.0_DP

    VelBulkCoefMin  =  0.0_DP
    TempBulkCoefMin =  0.0_DP
    QVapBulkCoefMin =  0.0_DP
    VelBulkCoefMax  =  1.0_DP
    TempBulkCoefMax =  1.0_DP
    QVapBulkCoefMax =  1.0_DP

    FlagFixFricTimeConstAtLB = .false.
    FricTimeConstAtLB        = 1.0e100_DP
    LowLatFricAtLB           = 1.0e100_DP
    FlagFixHeatFluxAtLB      = .false.
    HeatFluxAtLB             = 1.0e100_DP
    FlagFixMassFluxAtLB      = .false.
    MassFluxAtLB             = 1.0e100_DP

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, &                ! (in)
        & nml = surface_flux_bulk_nml, &  ! (out)
        & iostat = iostat_nml )        ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'BulkCoefMom', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'bulk coefficient for momentum', '1' )
    call HistoryAutoAddVariable( 'BulkCoefHeat', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'bulk coefficient for heat', '1' )
    call HistoryAutoAddVariable( 'SfcBulkRi', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'bulk Richardson number at the surface', '1' )

    call HistoryAutoAddVariable( 'TauX', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface stress(x)  ', 'N m-2' )
    call HistoryAutoAddVariable( 'TauY', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface stress(y)  ', 'N m-2' )
    call HistoryAutoAddVariable( 'Sens', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'sensible heat flux', 'W m-2' )
    call HistoryAutoAddVariable( 'SurfH2OVapFlux', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface H2O vapor flux  ', 'kg m-2 s-1' )
    call HistoryAutoAddVariable( 'Evap', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'latent heat flux  ', 'W m-2' )

    call HistoryAutoAddVariable( 'TauXB', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface stress(x)  ', 'N m-2' )
    call HistoryAutoAddVariable( 'TauYB', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface stress(y)  ', 'N m-2' )
    call HistoryAutoAddVariable( 'SensB', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'sensible heat flux', 'W m-2' )
    call HistoryAutoAddVariable( 'SurfH2OVapFluxB', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface H2O vapor flux  ', 'kg m-2 s-1' )
    call HistoryAutoAddVariable( 'EvapB', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'latent heat flux  ', 'W m-2' )

    call HistoryAutoAddVariable( 'TauXA', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface stress(x)  ', 'N m-2' )
    call HistoryAutoAddVariable( 'TauYA', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface stress(y)  ', 'N m-2' )
    call HistoryAutoAddVariable( 'SensA', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'sensible heat flux', 'W m-2' )
    call HistoryAutoAddVariable( 'SurfH2OVapFluxA', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface H2O vapor flux  ', 'kg m-2 s-1' )
    call HistoryAutoAddVariable( 'EvapA', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'latent heat flux  ', 'W m-2' )

    call HistoryAutoAddVariable( 'SurfH2OVapFluxU', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'surface H2O vapor flux  ', 'kg m-2 s-1' )
    call HistoryAutoAddVariable( 'EvapU', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'latent heat flux  ', 'W m-2' )

    call HistoryAutoAddVariable( 'MOLength', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'Monin-Obukhov length', 'm' )
    call HistoryAutoAddVariable( 'MOLengthInv', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'Monin-Obukhov length inverse', 'm-1' )
    call HistoryAutoAddVariable( 'BetaW', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'beta w_*', 'm s-1' )
    call HistoryAutoAddVariable( 'BLHeight', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'boundary layer height', 'm' )


    ! Initialization of modules used in this module
    !
    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    call SaturateInit


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )

    call MessageNotify( 'M', module_name, '  VelMinForRi       = %f', d = (/ VelMinForRi   /) )
    call MessageNotify( 'M', module_name, '  VelMinForVel      = %f', d = (/ VelMinForVel  /) )
    call MessageNotify( 'M', module_name, '  VelMinForTemp     = %f', d = (/ VelMinForTemp /) )
    call MessageNotify( 'M', module_name, '  VelMinForQVap     = %f', d = (/ VelMinForQVap /) )
    call MessageNotify( 'M', module_name, '  VelMaxForVel      = %f', d = (/ VelMaxForVel  /) )
    call MessageNotify( 'M', module_name, '  VelMaxForTemp     = %f', d = (/ VelMaxForTemp /) )
    call MessageNotify( 'M', module_name, '  VelMaxForQVap     = %f', d = (/ VelMaxForQVap /) )
    call MessageNotify( 'M', module_name, 'Bulk coefficients:' )
    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, '  QVapBulkCoefMin = %f', d = (/ QVapBulkCoefMin /) )
    call MessageNotify( 'M', module_name, '  VelBulkCoefMax  = %f', d = (/ VelBulkCoefMax  /) )
    call MessageNotify( 'M', module_name, '  TempBulkCoefMax = %f', d = (/ TempBulkCoefMax /) )
    call MessageNotify( 'M', module_name, '  QVapBulkCoefMax = %f', d = (/ QVapBulkCoefMax /) )
    call MessageNotify( 'M', module_name, 'FlagIncludeB94W          = %b', l = (/ FlagIncludeB94W /) )
    call MessageNotify( 'M', module_name, 'FlagFixFricTimeConstAtLB = %b', l = (/ FlagFixFricTimeConstAtLB /) )
    call MessageNotify( 'M', module_name, 'FricTimeConstAtLB        = %f', d = (/ FricTimeConstAtLB /) )
    call MessageNotify( 'M', module_name, 'LowLatFricAtLB           = %f', d = (/ LowLatFricAtLB /) )
    call MessageNotify( 'M', module_name, 'FlagFixHeatFluxAtLB      = %b', l = (/ FlagFixHeatFluxAtLB /) )
    call MessageNotify( 'M', module_name, 'HeatFluxAtLB             = %f', d = (/ HeatFluxAtLB /) )
    call MessageNotify( 'M', module_name, 'FlagFixMassFluxAtLB      = %b', l = (/ FlagFixMassFluxAtLB /) )
    call MessageNotify( 'M', module_name, 'MassFluxAtLB             = %f', d = (/ MassFluxAtLB /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    surface_flux_bulk_inited = .true.

  end subroutine SurfaceFluxInit

  !--------------------------------------------------------------------------------------

end module surface_flux_bulk