!= 地表面フラックス (バルク法)
!
!= Surface flux (Bulk method)
!
! Authors::   Masatsugu Odaka, Yoshiyuki O. Takahashi, Yasuhiro MORIKAWA, Yukiko YAMADA
! Version::   $Id: surfaceflux_bulk_core.f90,v 1.1 2012-10-01 06:25:03 odakker Exp $ 
! Tag Name::  $Name: arare5-20130116 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2012. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module surfaceflux_bulk_core
  !
  != 地表面フラックス
  !
  != Surface flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! Louis et al. (1982) の方法に基づいて地表面フラックスを計算. 
  !
  ! Surface fluxes are calculated by using the scheme by Louis et al. (1982). 
  !
  !== 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.
  !
  !== 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: ncmax   ! 化学種の数
                             ! Number og chemical scpecies

!!$  use gridset, only:
  use gridset_surfaceflux, 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 :: SurfaceFluxInit

  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  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
  real(DP), save:: Vel0 = 0.0d0  ! 下層での水平速度嵩上げ値
                                   ! 

  ! バルク係数
  ! 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: arare5-20130116 $' // &
    & '$Id: surfaceflux_bulk_core.f90,v 1.1 2012-10-01 06:25:03 odakker Exp $'
                              ! モジュールのバージョン
                              ! Module version


contains

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

  subroutine SurfaceFlux( &
    & xyz_U, xyz_V, xyz_Temp, xyr_Temp, xyr_VirTemp, xyzf_QMix,   & ! (in)
    & xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner, & ! (in)
    & xy_SurfTemp, xy_SurfHumidCoef,                              & ! (in)
    & xy_SurfRoughLength,                                         & ! (in)
    & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux,    & ! (inout)
    & xy_SurfVelTransCoef, xy_SurfTempTransCoef,                  & ! (out)
    & xy_SurfQVapTransCoef                                        & ! (out)
    & )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! Calculate radiation flux from temperature, specific humidity, and 
    ! air pressure. 
    !

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & 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
!!$      & LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation
      & TempSfc, &               
                              ! 地表面温度
                              ! Surface temperature

      & MolWtDry
                              ! 乾燥空気分子量 
                              ! Molecular waeight of dry air


    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
!!$    use saturate, only: xy_CalcQVapSat
    use composition, only: IdxCG,       & ! 凝結過程(気体)の配列添字 
                                          ! Index of vapor
      &                    IdxCC,       & ! 凝結過程(雲)の配列添字 
                                          ! Index of cloud
      &                    SpcWetID,    & !凝結成分の化学種ID
                                          ! ID number of moist condensation spices
      &                    CondNum,     & ! 雲の数
                                          ! Number of cloud component
      &                    MolWtWet,    & ! 凝結成分分子量
                                          ! Molecular weight of moist air
      &                    SpcWetSymbol   ! 凝結成分の化学種名
                                          ! Chemical formula of condensation spices

    ! 化学量計算に関する設定
    ! Settings for chemical calculation
    !
    use chemcalc, only: SvapPress ! 飽和蒸気圧
                                  ! Saturation vapor pressure

    ! 座標データ設定
    ! 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
    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):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ T_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    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_SurfRoughLength (0:imax-1, 1:jmax)
                              ! 地表粗度長. 
                              ! Surface rough length
    real(DP), intent(inout):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(inout):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(inout):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(inout):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    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

    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_SurfBulkRiNum (0:imax-1, 1:jmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $ number
    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_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

    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 )


    ! バルク $ R_i $ 数算出
    ! Calculate bulk $ R_i $
    !
    do j = 1, jmax
      do i = 0, imax-1
        xy_SurfVelAbs(i,j) = sqrt ( xyz_U(i,j,1)**2 + xyz_V(i,j,1)**2 ) + Vel0

        xy_SurfBulkRiNum(i,j) =                              &
          &   Grav / ( xy_SurfTemp(i,j) / xyr_Exner(i,j,0) ) &
          &   * ( xyz_Temp(i,j,1) / xyz_Exner(i,j,1)         &
          &     - xy_SurfTemp(i,j) / xyr_Exner(i,j,0) )      &
          &   / max( xy_SurfVelAbs(i,j), VelMinForRi )**2    &
          &   * ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) )
      end do
    end do


    ! バルク係数算出
    ! Calculate bulk coefficients
    !
    call BulkCoef( &
      & xy_SurfBulkRiNum,    & ! (in)
      & xy_SurfRoughLength,  & ! (in)
      & xy_SurfHeight,       & ! (in)
      & xyz_Height,          & ! (in)
      & xy_SurfVelBulkCoef,  & ! (out)
      & xy_SurfTempBulkCoef, & ! (out)
      & xy_SurfQVapBulkCoef  & ! (out)
      & )

    ! 輸送係数の計算
    ! 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_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) )


    ! 地表面フラックスの計算
    ! 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 )
      do n = 1, CondNum
        xyf_QMixFluxSurf(:,:,IdxCG(n)) = &
          & - xy_SurfHumidCoef * xy_SurfQVapTransCoef(:,:)           &
          & * ( xyzf_QMix(:,:,1,IdxCG(n))                            &
          &     - SvapPress( SpcWetID(IdxCC(n)), TempSfc )           &
          &             / ( xyr_Press(:,:,0) )                       &
          &             * (MolWtWet(IdxCG(n)) / MolWtDry)            &
          &  )
      end do
    else
      do n = 1, CondNum
        xyf_QMixFluxSurf(:,:,IdxCG(n)) = MassFluxAtLB
      end do
    end if
    !


!!$    xyf_QMixFluxSurf(:,:,1:IndexH2OVap-1)     = 0.0d0
!!$    xyf_QMixFluxSurf(:,:,IndexH2OVap+1:ncmax) = 0.0d0

    ! 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
    !
    xyr_MomFluxX(:,:,0) = xy_MomFluxXSurf(:,:)
    xyr_MomFluxY(:,:,0) = xy_MomFluxYSurf(:,:)
    xyr_HeatFlux(:,:,0) = xy_HeatFluxSurf(:,:)
    do n = 1, CondNum
      xyrf_QMixFlux(:,:,0,IdxCG(n)) = xyf_QMixFluxSurf(:,:,IdxCG(n))
    end do



    ! ヒストリデータ出力
    ! History data output
    !

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

  end subroutine SurfaceFlux

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

  subroutine BulkCoef( &
    & xy_SurfBulkRiNum, &       ! (in)
    & xy_SurfRoughLength,  & ! (in)
    & xy_SurfHeight,          & ! (in)
    & xyz_Height,             & ! (in)
    & xy_SurfVelBulkCoef,     & ! (out)
    & xy_SurfTempBulkCoef,    & ! (out)
    & xy_SurfQVapBulkCoef     & ! (out)
    & )
    !
    ! バルク係数を算出します.
    !
    ! Bulk coefficients are calculated.
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
!!$    use constants, only: &
    use constants0, only: &
      & FKarm                 ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xy_SurfBulkRiNum (0:imax-1, 1:jmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $ number
    real(DP), intent(in):: xy_SurfRoughLength (0:imax-1, 1:jmax)
                              ! 地表粗度長
                              ! Surface rough length
    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(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

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_SurfBulkCoefInNeutCond    (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

    else

      ! Parameterization by Louis et al. (1981)
      !

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

      if ( FlagUseOfBulkCoefInNeutralCond ) then

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

        xy_SurfVelBulkCoef  = xy_SurfBulkCoefInNeutCond
        xy_SurfTempBulkCoef = xy_SurfBulkCoefInNeutCond

        xy_SurfQVapBulkCoef = xy_SurfTempBulkCoef

      else

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

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

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

              xy_SurfVelBulkCoef(i,j) =                                        &
                &   xy_SurfBulkCoefInNeutCond(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_SurfBulkCoefInNeutCond(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_SurfBulkCoefInNeutCond(i,j)                                   &
                &   * (   1.0_DP                                                     &
                &       - 10.0_DP * xy_SurfBulkRiNum(i,j)                            &
                &         / (   1.0_DP                                               &
                &             + 75.0_DP * xy_SurfBulkCoefInNeutCond(i,j)             &
                &               * sqrt( - ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) ) &
                &                         / xy_SurfRoughLength(i,j)                  &
                &                         * xy_SurfBulkRiNum(i,j)                    &
                &                     )                                              &
                &           )                                                        &
                &     )

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

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

            end if
          end do
        end do

      end if

    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


  end subroutine BulkCoef

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

  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
    !

    ! 並列処理用変数
    ! Parallel processing variable 
    !       
    Use mpi_wrapper,only: myrank

    ! 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

    ! 宣言文 ; 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,                &
      !
      & VelMinForRi, VelMinForVel, VelMinForTemp, VelMinForQVap,      &
      & VelMaxForVel, VelMaxForTemp, VelMaxForQVap,                   &
      !
      & VelBulkCoefMin, TempBulkCoefMin, QVapBulkCoefMin,             &
      & VelBulkCoefMax, TempBulkCoefMax, QVapBulkCoefMax,             &
      !
      & FlagFixFricTimeConstAtLB, FricTimeConstAtLB, LowLatFricAtLB,  &
      & FlagFixHeatFluxAtLB,      HeatFluxAtLB,                       &
      & FlagFixMassFluxAtLB,      MassFluxAtLB, &
      !
      & Vel0
          !
          ! デフォルト値については初期化手続 "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

    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.0d100
    LowLatFricAtLB           = 1.0d100
    FlagFixHeatFluxAtLB      = .false.
    HeatFluxAtLB             = 1.0d100
    FlagFixMassFluxAtLB      = .false.
    MassFluxAtLB             = 1.0d100

    Vel0 = 0.0d0

    ! 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( 'TauX', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface stress(x)  ', 'N m-2' )
!!$    call HistoryAutoAddVariable( 'TauY', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface stress(y)  ', 'N m-2' )
!!$    call HistoryAutoAddVariable( 'Sens', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'sensible heat flux', 'W m-2' )
!!$    call HistoryAutoAddVariable( 'SurfH2OVapFlux', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface H2O vapor flux  ', 'kg m-2 s-1' )
!!$    call HistoryAutoAddVariable( 'Evap', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'latent heat flux  ', 'W m-2' )

!!$    call HistoryAutoAddVariable( 'TauXB', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface stress(x)  ', 'N m-2' )
!!$    call HistoryAutoAddVariable( 'TauYB', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface stress(y)  ', 'N m-2' )
!!$    call HistoryAutoAddVariable( 'SensB', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'sensible heat flux', 'W m-2' )
!!$    call HistoryAutoAddVariable( 'SurfH2OVapFluxB', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface H2O vapor flux  ', 'kg m-2 s-1' )
!!$    call HistoryAutoAddVariable( 'EvapB', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'latent heat flux  ', 'W m-2' )

!!$    call HistoryAutoAddVariable( 'TauXA', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface stress(x)  ', 'N m-2' )
!!$    call HistoryAutoAddVariable( 'TauYA', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface stress(y)  ', 'N m-2' )
!!$    call HistoryAutoAddVariable( 'SensA', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'sensible heat flux', 'W m-2' )
!!$    call HistoryAutoAddVariable( 'SurfH2OVapFluxA', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'surface H2O vapor flux  ', 'kg m-2 s-1' )
!!$    call HistoryAutoAddVariable( 'EvapA', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'latent heat flux  ', 'W m-2' )

    ! 印字 ; Print
    !
    if (myrank == 0) then 

    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, '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, "Vel0 = %f", d=(/Vel0/))
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
    end if

    surface_flux_bulk_inited = .true.

  end subroutine SurfaceFluxInit

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

end module surfaceflux_bulk_core
