Class surface_flux_bulk
In: surface_flux/surface_flux_bulk.f90

地表面フラックス

Surface flux

Note that Japanese and English are described in parallel.

地表面フラックスを計算.

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,

  1. Appl. Meteor., 30, 327-341, 1991.

Beljaars, A. C. M., The parameterization of surface fluxes in large-scale models under free convection,

  1. 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

Methods

Included Modules

gridset composition dc_types dc_message gtool_historyauto constants0 constants saturate axesset timeset dc_trace namelist_util dc_iounit dc_string

Public Instance methods

Subroutine :
BulkCoefMethod :character(*), intent(in)
: Method for calculating bulk coefficient
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)
xyz_VirTemp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T_v $ . 仮温度 (整数レベル). Virtual temperature (full level)
xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ T_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
xy_SurfVirTemp(0:imax-1, 1:jmax) :real(DP), intent(in)
: $ T_v $ . 仮温度 (惑星表面). Virtual temperature (surface)
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_SurfRoughLengthMom(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長. Surface rough length for momentum
xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長. Surface rough length for heat
xy_MomFluxX(0:imax-1, 1:jmax) :real(DP), intent(out)
: 惑星表面東西方向運動量フラックス. Eastward momentum flux at surface
xy_MomFluxY(0:imax-1, 1:jmax) :real(DP), intent(out)
: 惑星表面南北方向運動量フラックス. Northward momentum flux at surface
xy_HeatFlux(0:imax-1, 1:jmax) :real(DP), intent(out)
: 惑星表面熱フラックス. Heat flux at surface
xyf_QMixFlux(0:imax-1, 1:jmax, 1:ncmax) :real(DP), intent(out)
: 惑星表面比湿フラックス. Specific humidity flux at surface
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( BulkCoefMethod, xyz_U, xyz_V, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xy_SurfVirTemp, xyzf_QMix, xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfRoughLengthMom, xy_SurfRoughLengthHeat, xy_MomFluxX, xy_MomFluxY, xy_HeatFlux, xyf_QMixFlux, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! 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, CpDry, LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation

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

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimeN, 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(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

    ! 作業変数
    ! 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_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, xy_SurfRoughLengthMom , xy_SurfRoughLengthHeat, xy_SurfHeight, xyz_Height, xyz_U(:,:,1), xyz_V(:,:,1), xy_SurfVirTemp, xyz_VirTemp(:,:,1), xyr_Exner(:,:,0), xyz_Exner(:,:,1), xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef, xy_BetaW )


    ! 
    ! 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) )


    ! 地表面フラックスの計算
    ! 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.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
    !
    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 :

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
    !

    ! 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

    ! 宣言文 ; 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.0d100
    LowLatFricAtLB           = 1.0d100
    FlagFixHeatFluxAtLB      = .false.
    HeatFluxAtLB             = 1.0d100
    FlagFixMassFluxAtLB      = .false.
    MassFluxAtLB             = 1.0d100

    ! 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( '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' )

    ! 印字 ; 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
Subroutine :
xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 東西方向運動量フラックス. Eastward momentum flux
xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 南北方向運動量フラックス. Northward momentum flux
xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 熱フラックス. Heat flux
xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) :real(DP), intent(in)
: 比湿フラックス. Specific humidity flux
xy_SurfH2OVapFluxA(0:imax-1, 1:jmax) :real(DP), intent(in)
: 惑星表面水蒸気フラックス. Water vapor flux at the surface
xy_SurfLatentHeatFluxA(0:imax-1, 1:jmax) :real(DP), intent(in)
: 惑星表面潜熱フラックス. Latent heat flux at the surface
xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ DP{u}{t} $ . 東西風速時間変化率. Eastward wind tendency
xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ DP{v}{t} $ . 南北風速時間変化率. Northward wind tendency
xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ DP{T}{t} $ . 温度時間変化率. Temperature tendency
xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(in)
: $ DP{q}{t} $ . 比湿時間変化率. Specific humidity tendency
xy_SurfTemp(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表面温度. Surface temperature
xy_DSurfTempDt(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表面温度時間変化率. Surface temperature tendency
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half 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_SurfHumidCoef(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表湿潤度. Surface humidity coefficient
xy_SurfVelTransCoef(0:imax-1, 1:jmax) :real(DP), intent(in)
: 輸送係数:運動量. Diffusion coefficient: velocity
xy_SurfTempTransCoef(0:imax-1, 1:jmax) :real(DP), intent(in)
: 輸送係数:温度. Transfer coefficient: temperature
xy_SurfQVapTransCoef(0:imax-1, 1:jmax) :real(DP), intent(in)
: 輸送係数:水蒸気 Transfer coefficient: water vapor

フラックス (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.

[Source]

  subroutine SurfaceFluxOutput( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFluxA, xy_SurfLatentHeatFluxA, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xy_SurfTemp, xy_DSurfTempDt, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfHumidCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef )
    !
    ! フラックス (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, GasRDry, CpDry, LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation

    ! 飽和比湿の算出
    ! Evaluation of saturation specific humidity
    !
    use saturate, only: xy_CalcQVapSat, xy_CalcDQVapSatDTemp

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

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

    ! 宣言文 ; Declaration statements
    !
    implicit none
    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_SurfQVapSat (0:imax-1, 1:jmax)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity 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_SurfQVapSat       = xy_CalcQVapSat      ( xy_SurfTemp, xyr_Press(:,:,0) )
    xy_SurfDQVapSatDTemp = xy_CalcDQVapSatDTemp( xy_SurfTemp, xy_SurfQVapSat   )

    ! 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.0d0 * DelTime

        xyr_MomFluxYCor( i,j,0 ) = xyr_MomFluxY( i,j,0 ) - xy_SurfVelTransCoef( i,j ) * xyz_DVDt( i,j,1 ) * 2.0d0 * 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.0d0 * 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.0d0 * 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

Private Instance methods

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

Calculation of Psi_M and Psi_H

[Source]

  subroutine BH91CalcPsi( Gamma, xy_ZetaM, xy_ZetaH, xy_PsiM, xy_PsiH )
    !
    ! 
    !
    ! 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
          ParamXH = ( 1.0_DP - Gamma * ZetaH )**0.25

          ! 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 :
IDBulkCoefMethod :integer , intent(in)
xy_SurfRoughLengthMom(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長 Surface rough length for momentum
xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長 Surface rough length for heat
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
xy_U(0:imax-1, 1:jmax) :real(DP), intent(in)
: Eastward wind velocity at lowest level
xy_V(0:imax-1, 1:jmax) :real(DP), intent(in)
: Northward wind velocity at lowest level
xy_SurfVirTemp(0:imax-1, 1:jmax) :real(DP), intent(in)
: Surface virtual temperature
xy_VirTemp(0:imax-1, 1:jmax) :real(DP), intent(in)
: Virtual temperature at lowest layer
xy_SurfExner(0:imax-1, 1:jmax) :real(DP), intent(in)
: Exner function at the surface
xy_Exner(0:imax-1, 1:jmax) :real(DP), intent(in)
: Exner function at lowest layer
xy_SurfVelBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:運動量. Bulk coefficient: temperature
xy_SurfTempBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:温度. Bulk coefficient: temperature
xy_SurfQVapBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:比湿. Bulk coefficient: specific humidity
xy_BetaW(0:imax-1, 1:jmax) :real(DP), intent(out)
: "vertical velocity" (B94)

バルク係数を算出します.

Bulk coefficients are calculated.

[Source]

  subroutine BulkCoef( IDBulkCoefMethod, xy_SurfRoughLengthMom , xy_SurfRoughLengthHeat, xy_SurfHeight, xyz_Height, xy_U, xy_V, xy_SurfVirTemp, xy_VirTemp, xy_SurfExner, xy_Exner, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef, xy_BetaW )
    !
    ! バルク係数を算出します.
    !
    ! 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, 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)

    ! 作業変数
    ! 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, xy_SurfRoughLengthMom , xy_SurfRoughLengthHeat, xy_SurfHeight, xyz_Height, xy_SurfBulkCoefMomInNeutCond, xy_SurfBulkCoefHeatInNeutCond, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef )

        end if

        xy_BetaW = 0.0_DP

      case ( IDBulkCoefMethodBH91B94 )

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

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

      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 :
xy_SurfRoughLengthMom(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長 Surface rough length for momentum
xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長 Surface rough length for heat
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
xy_U(0:imax-1, 1:jmax) :real(DP), intent(in)
: Eastward wind velocity at lowest level
xy_V(0:imax-1, 1:jmax) :real(DP), intent(in)
: Northward wind velocity at lowest level
xy_SurfVirTemp(0:imax-1, 1:jmax) :real(DP), intent(in)
: Surface virtual temperature
xy_VirTemp(0:imax-1, 1:jmax) :real(DP), intent(in)
: Virtual temperature at lowest layer
xy_SurfExner(0:imax-1, 1:jmax) :real(DP), intent(in)
: Exner function at the surface
xy_Exner(0:imax-1, 1:jmax) :real(DP), intent(in)
: Exner function at lowest layer
xy_SurfVelBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:運動量. Bulk coefficient: temperature
xy_SurfTempBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:温度. Bulk coefficient: temperature
xy_SurfQVapBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:比湿. Bulk coefficient: specific humidity
xy_BetaW(0:imax-1, 1:jmax) :real(DP), intent(out)
: "vertical velocity" (B94)
xy_SurfBulkRiNum(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク $ R_i $ 数. Bulk $ R_i $ number

バルク係数を算出します.

Bulk coefficients are calculated.

[Source]

  subroutine BulkCoefBH91B94( xy_SurfRoughLengthMom , xy_SurfRoughLengthHeat, xy_SurfHeight, xyz_Height, xy_U, xy_V, xy_SurfVirTemp, xy_VirTemp, xy_SurfExner, xy_Exner, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef, xy_BetaW, xy_SurfBulkRiNum )
    !
    ! バルク係数を算出します.
    !
    ! 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, 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

    ! 作業変数
    ! 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)
    logical :: a_FlagReCalcLocal (1)
    logical :: a_FlagReCalcGlobal(1)
    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.0d10
    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, xy_ZetaM, xy_ZetaH, xy_PsiM0, xy_PsiH0 )
        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, xy_ZetaM, xy_ZetaH, xy_PsiM1, xy_PsiH1 )
      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.0d-10
            xy_MOLength(i,j) = ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) ) / ( FKarm * SurfBulkRiNum ) * xy_SurfVelBulkCoef(i,j)**1.5 / 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.0d-10 )
            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


    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 :
xy_SurfBulkRiNum(0:imax-1, 1:jmax) :real(DP), intent(in)
: バルク $ R_i $ 数. Bulk $ R_i $ number
xy_SurfRoughLengthMom(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長 Surface rough length for momentum
xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) :real(DP), intent(in)
: 地表粗度長 Surface rough length for heat
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
xy_SurfBulkCoefMomInNeutCond(0:imax-1, 1:jmax) :real(DP), intent(in)
xy_SurfBulkCoefHeatInNeutCond(0:imax-1, 1:jmax) :real(DP), intent(in)
xy_SurfVelBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:運動量. Bulk coefficient: temperature
xy_SurfTempBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:温度. Bulk coefficient: temperature
xy_SurfQVapBulkCoef(0:imax-1, 1:jmax) :real(DP), intent(out)
: バルク係数:比湿. Bulk coefficient: specific humidity

バルク係数を算出します.

Bulk coefficients are calculated.

[Source]

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

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

    ! 宣言文 ; 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

    ! 作業変数
    ! Work variables
    !

    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
      end do
    end do


  end subroutine BulkCoefL82
ConstBulkCoef
Variable :
ConstBulkCoef :real(DP), save
: バルク係数一定値. Steady value of bulk coefficient
FlagConstBulkCoef
Variable :
FlagConstBulkCoef :logical, save
: Flag for using constant bulk coefficient
FlagFixFricTimeConstAtLB
Variable :
FlagFixFricTimeConstAtLB :logical , save
FlagFixHeatFluxAtLB
Variable :
FlagFixHeatFluxAtLB :logical , save
FlagFixMassFluxAtLB
Variable :
FlagFixMassFluxAtLB :logical , save
FlagIncludeB94W
Variable :
FlagIncludeB94W :logical
FlagUseOfBulkCoefInNeutralCond
Variable :
FlagUseOfBulkCoefInNeutralCond :logical, save
: Flag for using bulk coefficient in neutral condition
FricTimeConstAtLB
Variable :
FricTimeConstAtLB :real(DP), save
: 下部境界摩擦の時定数 (s). Time constant of surface friction (s).
HeatFluxAtLB
Variable :
HeatFluxAtLB :real(DP), save
: 下部境界での熱フラックス (W m-2). Heat flux at the lower boundary (W m-2).
IDBulkCoefMethodBH91B94
Constant :
IDBulkCoefMethodBH91B94 = 2 :integer, parameter
IDBulkCoefMethodL82
Constant :
IDBulkCoefMethodL82 = 1 :integer, parameter
LowLatFricAtLB
Variable :
LowLatFricAtLB :real(DP), save
: 下部境界摩擦が働く最低緯度 (degree). Lowest latitude where the friction is applied (degree)
MassFluxAtLB
Variable :
MassFluxAtLB :real(DP), save
: 下部境界での質量フラックス (W m-2). 実際にはゼロに固定するために使う程度にしか使えないだろう. Mass flux at the lower boundary (kg m-2 s-1).
QVapBulkCoefMax
Variable :
QVapBulkCoefMax :real(DP), save
: $ q $ バルク係数最大値. Maximum value of $ q $ bulk coefficient
QVapBulkCoefMin
Variable :
QVapBulkCoefMin :real(DP), save
: $ q $ バルク係数最小値. Minimum value of $ q $ bulk coefficient
TempBulkCoefMax
Variable :
TempBulkCoefMax :real(DP), save
: $ T $ バルク係数最大値. Maximum value of $ T $ bulk coefficient
TempBulkCoefMin
Variable :
TempBulkCoefMin :real(DP), save
: $ T $ バルク係数最小値. Minimum value of $ T $ bulk coefficient
VelBulkCoefMax
Variable :
VelBulkCoefMax :real(DP), save
: $ u $ バルク係数最大値. Maximum value of $ u $ bulk coefficient
VelBulkCoefMin
Variable :
VelBulkCoefMin :real(DP), save
: $ u $ バルク係数最小値. Minimum value of $ u $ bulk coefficient
VelMaxForQVap
Variable :
VelMaxForQVap :real(DP), save
: 水蒸気用風最大値. Maximum value of velocity for vapor
VelMaxForTemp
Variable :
VelMaxForTemp :real(DP), save
: 熱用風最大値. Maximum value of velocity for thermal
VelMaxForVel
Variable :
VelMaxForVel :real(DP), save
: 運動量用風最大値. Maximum value of velocity for momentum
VelMinForQVap
Variable :
VelMinForQVap :real(DP), save
: 水蒸気用風最小値. Minimum value of velocity for vapor
VelMinForRi
Variable :
VelMinForRi :real(DP), save
: $ R_i $ 数用風最小値. Minimum value of velocity for $ R_i $ number
VelMinForTemp
Variable :
VelMinForTemp :real(DP), save
: 熱用風最小値. Minimum value of velocity for thermal
VelMinForVel
Variable :
VelMinForVel :real(DP), save
: 運動量用風最小値. Minimum value of velocity for momentum
module_name
Constant :
module_name = ‘surface_flux_bulk :character(*), parameter
: モジュールの名称. Module name
surface_flux_bulk_inited
Variable :
surface_flux_bulk_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag
version
Constant :
version = ’$Name: dcpam5-20130921 $’ // ’$Id: surface_flux_bulk.f90,v 1.22 2013-09-21 14:38:43 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version