Class surfaceflux_bulk_core
In: ../src/surface_flux/surfaceflux_bulk_core.f90

地表面フラックス

Surface flux

Note that Japanese and English are described in parallel.

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

Methods

Included Modules

gridset gridset_surfaceflux dc_types dc_message constants0 constants composition chemcalc dc_trace mpi_wrapper namelist_util dc_iounit dc_string

Public Instance methods

Subroutine :
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ u $ . 東西風速. Eastward wind
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ v $ . 南北風速. Northward wind
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T $ . 温度 (整数レベル). Temperature (full level)
xyr_Temp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ T $ . 温度 (半整数レベル). Temperature (half level)
xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ T_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(in)
: $ q $ . 比湿. Specific humidity
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ p_s $ . 地表面気圧 (半整数レベル). Surface pressure (half level)
xy_SurfHeight(0:imax-1,1:jmax) :real(DP), intent(in)
: $ z_s $ . 地表面高度. Surface height.
xyz_Height(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full level)
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
xy_SurfTemp(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表面温度. Surface temperature
xy_SurfHumidCoef(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表湿潤度. Surface humidity coefficient
xy_SurfRoughLength(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長. Surface rough length
xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(inout)
: 東西方向運動量フラックス. Eastward momentum flux
xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(inout)
: 南北方向運動量フラックス. Northward momentum flux
xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(inout)
: 熱フラックス. Heat flux
xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) :real(DP), intent(inout)
: 比湿フラックス. Specific humidity flux
xy_SurfVelTransCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: 輸送係数:運動量. Diffusion coefficient: velocity
xy_SurfTempTransCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: 輸送係数:温度. Transfer coefficient: temperature
xy_SurfQVapTransCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: 輸送係数:水蒸気 Transfer coefficient: water vapor

温度, 比湿, 気圧から, 放射フラックスを計算します.

Calculate radiation flux from temperature, specific humidity, and air pressure.

[Source]

  subroutine SurfaceFlux( xyz_U, xyz_V, xyz_Temp, xyr_Temp, xyr_VirTemp, xyzf_QMix, xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfRoughLength, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! 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, GasRDry, CpDry, TempSfc, MolWtDry
                              ! 乾燥空気分子量 
                              ! Molecular waeight of dry air


    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
!!$    use saturate, only: xy_CalcQVapSat
    use composition, only: IdxCG, IdxCC, SpcWetID, CondNum, MolWtWet, 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, xy_SurfRoughLength, xy_SurfHeight, xyz_Height, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef )

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

surface_flux_bulk モジュールの初期化を行います. NAMELIST#surface_flux_bulk_nml の読み込みはこの手続きで行われます.

"surface_flux_bulk" module is initialized. "NAMELIST#surface_flux_bulk_nml" is loaded in this procedure.

This procedure input/output NAMELIST#surface_flux_bulk_nml .

[Source]

  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, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = surface_flux_bulk_nml, 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