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

BH91CalcPsi   BulkCoef   BulkCoefBH91B94   BulkCoefL82   ConstBulkCoef   FlagConstBulkCoef   FlagFixFricTimeConstAtLB   FlagFixHeatFluxAtLB   FlagFixMassFluxAtLB   FlagIncludeB94W   FlagUseOfBulkCoefInNeutralCond   FricTimeConstAtLB   HeatFluxAtLB   IDBulkCoefMethodBH91B94   IDBulkCoefMethodL82   LowLatFricAtLB   MassFluxAtLB   QVapBulkCoefMax   QVapBulkCoefMin   SurfaceFlux   SurfaceFluxInit   SurfaceFluxOutput   TempBulkCoefMax   TempBulkCoefMin   VelBulkCoefMax   VelBulkCoefMin   VelMaxForQVap   VelMaxForTemp   VelMaxForVel   VelMinForQVap   VelMinForRi   VelMinForTemp   VelMinForVel   module_name   surface_flux_bulk_inited   version  

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

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 )
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
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure

    ! 飽和比湿の算出
    ! 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
SurfaceFluxInit( )
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

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


    ! 宣言文 ; Declaration statements
    !
    implicit none

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

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

    ! 実行文 ; Executable statement
    !

    if ( surface_flux_bulk_inited ) return


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

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

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

    FlagFixFricTimeConstAtLB = .false.
    FricTimeConstAtLB        = 1.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' )


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


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

    call MessageNotify( 'M', module_name, '  VelMinForRi       = %f', d = (/ VelMinForRi   /) )
    call MessageNotify( 'M', module_name, '  VelMinForVel      = %f', d = (/ VelMinForVel  /) )
    call MessageNotify( 'M', module_name, '  VelMinForTemp     = %f', d = (/ VelMinForTemp /) )
    call MessageNotify( 'M', module_name, '  VelMinForQVap     = %f', d = (/ VelMinForQVap /) )
    call MessageNotify( 'M', module_name, '  VelMaxForVel      = %f', d = (/ VelMaxForVel  /) )
    call MessageNotify( 'M', module_name, '  VelMaxForTemp     = %f', d = (/ VelMaxForTemp /) )
    call MessageNotify( 'M', module_name, '  VelMaxForQVap     = %f', d = (/ VelMaxForQVap /) )
    call MessageNotify( 'M', module_name, 'Bulk coefficients:' )
    call MessageNotify( 'M', module_name, '  FlagConstBulkCoef              = %b', l = (/ FlagConstBulkCoef /) )
    call MessageNotify( 'M', module_name, '  FlagUseOfBulkCoefInNeutralCond = %b', l = (/ FlagUseOfBulkCoefInNeutralCond /) )
    call MessageNotify( 'M', module_name, '  ConstBulkCoef   = %f', d = (/ ConstBulkCoef   /) )
    call MessageNotify( 'M', module_name, '  VelBulkCoefMin  = %f', d = (/ VelBulkCoefMin  /) )
    call MessageNotify( 'M', module_name, '  TempBulkCoefMin = %f', d = (/ TempBulkCoefMin /) )
    call MessageNotify( 'M', module_name, '  QVapBulkCoefMin = %f', d = (/ QVapBulkCoefMin /) )
    call MessageNotify( 'M', module_name, '  VelBulkCoefMax  = %f', d = (/ VelBulkCoefMax  /) )
    call MessageNotify( 'M', module_name, '  TempBulkCoefMax = %f', d = (/ TempBulkCoefMax /) )
    call MessageNotify( 'M', module_name, '  QVapBulkCoefMax = %f', d = (/ QVapBulkCoefMax /) )
    call MessageNotify( 'M', module_name, 'FlagIncludeB94W          = %b', l = (/ FlagIncludeB94W /) )
    call MessageNotify( 'M', module_name, 'FlagFixFricTimeConstAtLB = %b', l = (/ FlagFixFricTimeConstAtLB /) )
    call MessageNotify( 'M', module_name, 'FricTimeConstAtLB        = %f', d = (/ FricTimeConstAtLB /) )
    call MessageNotify( 'M', module_name, 'LowLatFricAtLB           = %f', d = (/ LowLatFricAtLB /) )
    call MessageNotify( 'M', module_name, 'FlagFixHeatFluxAtLB      = %b', l = (/ FlagFixHeatFluxAtLB /) )
    call MessageNotify( 'M', module_name, 'HeatFluxAtLB             = %f', d = (/ HeatFluxAtLB /) )
    call MessageNotify( 'M', module_name, 'FlagFixMassFluxAtLB      = %b', l = (/ FlagFixMassFluxAtLB /) )
    call MessageNotify( 'M', module_name, 'MassFluxAtLB             = %f', d = (/ MassFluxAtLB /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    surface_flux_bulk_inited = .true.

  end subroutine SurfaceFluxInit
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 )
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

BH91CalcPsi( Gamma, xy_ZetaM, xy_ZetaH, xy_PsiM, xy_PsiH )
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
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 )
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
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 )
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
BulkCoefL82( xy_SurfBulkRiNum, xy_SurfRoughLengthMom, xy_SurfRoughLengthHeat, xy_SurfHeight, xyz_Height, xy_SurfBulkCoefMomInNeutCond, xy_SurfBulkCoefHeatInNeutCond, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef )
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-20140218 $’ // ’$Id: surface_flux_bulk.f90,v 1.24 2014-02-14 08:12:20 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version