| Class | surface_flux_bulk |
| In: |
surface_flux/surface_flux_bulk.f90
|
Note that Japanese and English are described in parallel.
地表面フラックスを計算.
Surface fluxes are calculated.
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,
Beljaars, A. C. M., The parameterization of surface fluxes in large-scale models under free convection,
| SurfaceFlux : | 地表面フラックスの計算 |
| SurfaceFluxOutput : | 地表面フラックスの出力 |
| ———— : | ———— |
| SurfaceFlux : | Calculate surface fluxes |
| SurfaceFluxOutput : | Output surface fluxes |
| Subroutine : | |||
| BulkCoefMethod : | character(*), intent(in)
| ||
| xyz_U(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_V(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Temp(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_VirTemp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfVirTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfHeight(0:imax-1,1:jmax) : | real(DP), intent(in)
| ||
| xyz_Height(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHumidCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfRoughLengthMom(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SnowFrac(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_MomFluxX(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_MomFluxY(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_HeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xyf_QMixFlux(0:imax-1, 1:jmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xy_SurfVelTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfQVapTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfMOLength(0:imax-1, 1:jmax) : | real(DP), intent(out) |
温度, 比湿, 気圧から, 放射フラックスを計算します.
Calculate radiation flux from temperature, specific humidity, and air pressure.
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_SnowFrac, xy_MomFluxX, xy_MomFluxY, xy_HeatFlux, xyf_QMixFlux, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xy_SurfMOLength )
!
! 温度, 比湿, 気圧から, 放射フラックスを計算します.
!
! 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_CalcQVapSatOnLiq, xy_CalcQVapSatOnSol
! 座標データ設定
! 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(in):: xy_SnowFrac(0:imax-1, 1:jmax)
!
! Snow fraction
real(DP), intent(out):: xy_MomFluxX (0:imax-1, 1:jmax)
! 惑星表面東西方向運動量フラックス.
! Eastward momentum flux at surface
real(DP), intent(out):: xy_MomFluxY (0:imax-1, 1:jmax)
! 惑星表面南北方向運動量フラックス.
! Northward momentum flux at surface
real(DP), intent(out):: xy_HeatFlux (0:imax-1, 1:jmax)
! 惑星表面熱フラックス.
! Heat flux at surface
real(DP), intent(out):: xyf_QMixFlux(0:imax-1, 1:jmax, 1:ncmax)
! 惑星表面比湿フラックス.
! Specific humidity flux at surface
real(DP), intent(out):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
! 輸送係数:運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数:温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
! 輸送係数:水蒸気
! Transfer coefficient: water vapor
real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)
! 作業変数
! Work variables
!
real(DP):: xy_SurfTempBulkCoef (0:imax-1, 1:jmax)
! バルク係数:温度.
! Bulk coefficient: temperature
real(DP):: xy_SurfQVapBulkCoef (0:imax-1, 1:jmax)
! バルク係数:比湿.
! Bulk coefficient: specific humidity
real(DP):: xy_SurfVelBulkCoef (0:imax-1, 1:jmax)
! バルク係数:運動量.
! Bulk coefficient: temperature
real(DP):: xy_SurfVelAbs (0:imax-1, 1:jmax)
! 風速絶対値.
! Absolute velocity
real(DP):: xy_SurfQVapSatOnLiq (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_MomFluxXSurf (0:imax-1, 1:jmax)
! 地表面の東西方向運動量フラックス.
! Eastward momentum flux on surface
real(DP):: xy_MomFluxYSurf (0:imax-1, 1:jmax)
! 地表面の南北方向運動量フラックス.
! Northward momentum flux on surface
real(DP):: xy_HeatFluxSurf (0:imax-1, 1:jmax)
! 地表面の熱フラックス.
! Heat flux on surface
real(DP):: xyf_QMixFluxSurf(0:imax-1, 1:jmax, 1:ncmax)
! 地表面の質量フラックス.
! Mass flux of constituents on surface
real(DP):: xy_BetaW (0:imax-1, 1:jmax)
!
! "vertical velocity" (B94)
integer :: IDBulkCoefMethod
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. surface_flux_bulk_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! Check method for calculating bulk coefficient
!
if ( BulkCoefMethod == 'L82' ) then
IDBulkCoefMethod = IDBulkCoefMethodL82
else if ( BulkCoefMethod == 'BH91B94' ) then
IDBulkCoefMethod = IDBulkCoefMethodBH91B94
else
call MessageNotify( 'E', module_name, 'BulkCoefMethod of %c is inappropriate.', c1 = trim( BulkCoefMethod ) )
end if
! バルク係数算出
! Calculate bulk coefficients
!
call BulkCoef( IDBulkCoefMethod, 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, xy_SurfMOLength )
!
! Calculation of wind speed
!
xy_SurfVelAbs = sqrt ( xyz_U(:,:,1)**2 + xyz_V(:,:,1)**2 + xy_BetaW**2 )
!!$ xy_SurfVelAbs = sqrt ( xyz_U(:,:,1)**2 + xyz_V(:,:,1)**2 )
! 輸送係数の計算
! Calculate transfer coefficient
!
if ( .not. FlagFixFricTimeConstAtLB ) then
do i = 0, imax-1
do j = 1, jmax
!!$ xy_SurfVelTransCoef(i,j) = &
!!$ & xy_SurfVelBulkCoef(i,j) &
!!$ & * xyr_Press(i,j,0) / ( GasRDry * xyr_Temp(i,j,0) ) &
!!$ & * min( max( xy_SurfVelAbs(i,j), VelMinForVel ), VelMaxForVel )
xy_SurfVelTransCoef(i,j) = xy_SurfVelBulkCoef(i,j) * xyr_Press(i,j,0) / ( GasRDry * xyr_VirTemp(i,j,0) ) * min( max( xy_SurfVelAbs(i,j), VelMinForVel ), VelMaxForVel )
end do
end do
else
do j = 1, jmax
if ( abs( y_Lat(j) ) >= LowLatFricAtLB * PI / 180.0_DP ) then
xy_SurfVelTransCoef(:,j) = 1.0_DP / FricTimeConstAtLB
else
xy_SurfVelTransCoef(:,j) = 0.0_DP
end if
end do
end if
if ( .not. FlagFixHeatFluxAtLB ) then
do i = 0, imax-1
do j = 1, jmax
!!$ xy_SurfTempTransCoef(i,j) = &
!!$ & xy_SurfTempBulkCoef(i,j) &
!!$ & * xyr_Press(i,j,0) / ( GasRDry * xyr_Temp(i,j,0) ) &
!!$ & * min( max( xy_SurfVelAbs(i,j), VelMinForTemp ), VelMaxForTemp )
xy_SurfTempTransCoef(i,j) = xy_SurfTempBulkCoef(i,j) * xyr_Press(i,j,0) / ( GasRDry * xyr_VirTemp(i,j,0) ) * min( max( xy_SurfVelAbs(i,j), VelMinForTemp ), VelMaxForTemp )
end do
end do
else
! Set meaningless value.
xy_SurfTempTransCoef = 0.0_DP
end if
if ( .not. FlagFixMassFluxAtLB ) then
do i = 0, imax-1
do j = 1, jmax
!!$ xy_SurfQVapTransCoef(i,j) = &
!!$ & xy_SurfQVapBulkCoef(i,j) &
!!$ & * xyr_Press(i,j,0) / ( GasRDry * xyr_Temp(i,j,0) ) &
!!$ & * min( max( xy_SurfVelAbs(i,j), VelMinForQVap ), VelMaxForQVap )
xy_SurfQVapTransCoef(i,j) = xy_SurfQVapBulkCoef(i,j) * xyr_Press(i,j,0) / ( GasRDry * xyr_VirTemp(i,j,0) ) * min( max( xy_SurfVelAbs(i,j), VelMinForQVap ), VelMaxForQVap )
end do
end do
else
! Set meaningless value.
xy_SurfQVapTransCoef = 0.0_DP
end if
! 飽和比湿の計算
! Calculate saturated specific humidity
!
!!$ xy_SurfQVapSat = xy_CalcQVapSat( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSatOnLiq = xy_CalcQVapSatOnLiq( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSatOnSol = xy_CalcQVapSatOnSol( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSat = ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq + xy_SnowFrac * xy_SurfQVapSatOnSol
! 地表面フラックスの計算
! Calculate fluxes on flux
!
! Momentum
!
xy_MomFluxXSurf = - xy_SurfVelTransCoef * xyz_U(:,:,1)
xy_MomFluxYSurf = - xy_SurfVelTransCoef * xyz_V(:,:,1)
! Heat
!
if ( .not. FlagFixHeatFluxAtLB ) then
xy_HeatFluxSurf = - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef * ( xyz_Temp(:,:,1) / xyz_Exner(:,:,1) - xy_SurfTemp / xyr_Exner(:,:,0) )
else
xy_HeatFluxSurf = HeatFluxAtLB
end if
! Mass
!
if ( .not. FlagFixMassFluxAtLB ) then
xyf_QMixFluxSurf(:,:,IndexH2OVap) = - xy_SurfHumidCoef * xy_SurfQVapTransCoef(:,:) * ( xyzf_QMix(:,:,1,IndexH2OVap) - xy_SurfQVapSat )
else
xyf_QMixFluxSurf(:,:,IndexH2OVap) = MassFluxAtLB
end if
!
xyf_QMixFluxSurf(:,:,1:IndexH2OVap-1) = 0.0_DP
xyf_QMixFluxSurf(:,:,IndexH2OVap+1:ncmax) = 0.0_DP
! Surface flux of constituents except for water vapor is zero.
!!$ write( 6, * ) "MEMO: Surface flux of constituents except for water vapor is zero. (YOT, 2009/08/14)"
! フラックスの計算
! Calculate fluxes
!
xy_MomFluxX = xy_MomFluxXSurf
xy_MomFluxY = xy_MomFluxYSurf
xy_HeatFlux = xy_HeatFluxSurf
do n = 1, ncmax
xyf_QMixFlux(:,:,n) = xyf_QMixFluxSurf(:,:,n)
end do
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'BulkCoefMom' , xy_SurfVelBulkCoef )
call HistoryAutoPut( TimeN, 'BulkCoefHeat', xy_SurfTempBulkCoef )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine SurfaceFlux
| Subroutine : |
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 .
subroutine SurfaceFluxInit
!
! surface_flux_bulk モジュールの初期化を行います.
! NAMELIST#surface_flux_bulk_nml の読み込みはこの手続きで行われます.
!
! "surface_flux_bulk" module is initialized.
! "NAMELIST#surface_flux_bulk_nml" is loaded in this procedure.
!
! モジュール引用 ; USE statements
!
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output
! 文字列操作
! Character handling
!
use dc_string, only: StoA
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoAddVariable
! 座標データ設定
! Axes data settings
!
use axesset, only: AxnameX, AxnameY, AxnameZ, AxnameR, AxnameT
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
use saturate, only: SaturateInit
! 宣言文 ; Declaration statements
!
implicit none
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
namelist /surface_flux_bulk_nml/ FlagConstBulkCoef, FlagUseOfBulkCoefInNeutralCond, ConstBulkCoef, FlagIncludeB94W, VelMinForRi, VelMinForVel, VelMinForTemp, VelMinForQVap, VelMaxForVel, VelMaxForTemp, VelMaxForQVap, VelBulkCoefMin, TempBulkCoefMin, QVapBulkCoefMin, VelBulkCoefMax, TempBulkCoefMax, QVapBulkCoefMax, FlagFixFricTimeConstAtLB, FricTimeConstAtLB, LowLatFricAtLB, FlagFixHeatFluxAtLB, HeatFluxAtLB, FlagFixMassFluxAtLB, MassFluxAtLB
!
! デフォルト値については初期化手続 "surface_flux_bulk#SurfaceFluxInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "surface_flux_bulk#SurfaceFluxInit" for the default values.
!
! 実行文 ; Executable statement
!
if ( surface_flux_bulk_inited ) return
! デフォルト値の設定
! Default values settings
!
FlagConstBulkCoef = .false.
FlagUseOfBulkCoefInNeutralCond = .false.
ConstBulkCoef = 0.0_DP
!!$ FlagIncludeB94W = .false.
FlagIncludeB94W = .true.
VelMinForRi = 0.01_DP
VelMinForVel = 0.01_DP
VelMinForTemp = 0.01_DP
VelMinForQVap = 0.01_DP
VelMaxForVel = 1000.0_DP
VelMaxForTemp = 1000.0_DP
VelMaxForQVap = 1000.0_DP
VelBulkCoefMin = 0.0_DP
TempBulkCoefMin = 0.0_DP
QVapBulkCoefMin = 0.0_DP
VelBulkCoefMax = 1.0_DP
TempBulkCoefMax = 1.0_DP
QVapBulkCoefMax = 1.0_DP
FlagFixFricTimeConstAtLB = .false.
FricTimeConstAtLB = 1.0e100_DP
LowLatFricAtLB = 1.0e100_DP
FlagFixHeatFluxAtLB = .false.
HeatFluxAtLB = 1.0e100_DP
FlagFixMassFluxAtLB = .false.
MassFluxAtLB = 1.0e100_DP
! NAMELIST の読み込み
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, 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
| Subroutine : | |||
| xy_SnowFrac(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_SurfH2OVapFluxA(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfLatentHeatFluxA(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfHumidCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfVelTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfQVapTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
|
フラックス (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux). について, その他の引数を用いて補正し, 出力を行う.
Fluxes (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux) are corrected by using other arguments, and the corrected values are output.
subroutine SurfaceFluxOutput( xy_SnowFrac, 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_CalcQVapSatOnLiq, xy_CalcQVapSatOnSol, xy_CalcDQVapSatDTempOnLiq, xy_CalcDQVapSatDTempOnSol
! 時刻管理
! 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):: xy_SnowFrac(0:imax-1, 1:jmax)
!
! Snow fraction
real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 比湿フラックス.
! Specific humidity flux
real(DP), intent(in):: xy_SurfH2OVapFluxA (0:imax-1, 1:jmax)
! 惑星表面水蒸気フラックス.
! Water vapor flux at the surface
real(DP), intent(in):: xy_SurfLatentHeatFluxA(0:imax-1, 1:jmax)
! 惑星表面潜熱フラックス.
! Latent heat flux at the surface
real(DP), intent(in):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{u}{t} $ . 東西風速時間変化率.
! Eastward wind tendency
real(DP), intent(in):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{v}{t} $ . 南北風速時間変化率.
! Northward wind tendency
real(DP), intent(in):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度時間変化率.
! Temperature tendency
real(DP), intent(in):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 比湿時間変化率.
! Specific humidity tendency
real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度時間変化率.
! Surface temperature tendency
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
! 地表湿潤度.
! Surface humidity coefficient
real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
! 輸送係数:運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数:温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
! 輸送係数:水蒸気
! Transfer coefficient: water vapor
! 出力のための作業変数
! Work variables for output
!
real(DP):: xyr_MomFluxXCor (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP):: xyr_MomFluxYCor (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP):: xyr_HeatFluxCor (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP):: xyrf_QMixFluxCor(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 比湿フラックス.
! Specific humidity flux
real(DP):: xyr_LatentHeatFluxCor(0:imax-1, 1:jmax, 0:kmax)
! 表面潜熱フラックス.
! Latent heat flux
real(DP):: xy_SurfQVapSatOnLiq (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
! 作業変数
! Work variables
!
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. surface_flux_bulk_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! 飽和比湿の計算
! Calculate saturated specific humidity
!
xy_SurfQVapSatOnLiq = xy_CalcQVapSatOnLiq( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSatOnSol = xy_CalcQVapSatOnSol( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSat = ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq + xy_SnowFrac * xy_SurfQVapSatOnSol
xy_SurfDQVapSatDTempOnLiq = xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
xy_SurfDQVapSatDTempOnSol = xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
xy_SurfDQVapSatDTemp = ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq + xy_SnowFrac * xy_SurfDQVapSatDTempOnSol
! Output of fluxes at t
!
! 風速, 温度, 比湿フラックス補正
! Correct fluxes of wind, temperature, specific humidity
!
do j = 1, jmax
do i = 0, imax-1
xyr_MomFluxXCor( i,j,0 ) = xyr_MomFluxX( i,j,0 ) - xy_SurfVelTransCoef( i,j ) * xyz_DUDt( i,j,1 ) * DelTime
xyr_MomFluxYCor( i,j,0 ) = xyr_MomFluxY( i,j,0 ) - xy_SurfVelTransCoef( i,j ) * xyz_DVDt( i,j,1 ) * DelTime
xyr_HeatFluxCor( i,j,0 ) = xyr_HeatFlux( i,j,0 ) - CpDry * xyr_Exner( i,j,0 ) * xy_SurfTempTransCoef( i,j ) * ( xyz_DTempDt( i,j,1 ) / xyz_Exner( i,j,1 ) - xy_DSurfTempDt( i,j ) / xyr_Exner( i,j,0 ) ) * DelTime
end do
end do
n = IndexH2OVap
do j = 1, jmax
do i = 0, imax-1
xyrf_QMixFluxCor( i,j,0,n ) = xyrf_QMixFlux( i,j,0,n ) - xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j ) * ( xyzf_DQMixDt( i,j,1,n ) - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * DelTime
end do
end do
do n = 1, IndexH2OVap-1
xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
end do
do n = IndexH2OVap+1, ncmax
xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
end do
n = IndexH2OVap
do j = 1, jmax
do i = 0, imax-1
!!$ xyr_LatentHeatFluxCor( i,j,0 ) = xy_SurfLatentHeatFlux( i,j ) &
!!$ & - LatentHeat * xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j ) &
!!$ & * ( xyzf_DQMixDt( i,j,1,n ) &
!!$ & - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * DelTime
xyr_LatentHeatFluxCor( i,j,0 ) = LatentHeat * xyrf_QMixFluxCor( i,j,0,n )
end do
end do
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'TauX' , xyr_MomFluxXCor (:,:,0) )
call HistoryAutoPut( TimeN, 'TauY' , xyr_MomFluxYCor (:,:,0) )
call HistoryAutoPut( TimeN, 'Sens' , xyr_HeatFluxCor (:,:,0) )
call HistoryAutoPut( TimeN, 'SurfH2OVapFlux', xyrf_QMixFluxCor(:,:,0,IndexH2OVap) )
call HistoryAutoPut( TimeN, 'Evap' , xyr_LatentHeatFluxCor(:,:,0) )
! Output of fluxes at t - \Delta t
!
! 風速, 温度, 比湿フラックス補正
! Correct fluxes of wind, temperature, specific humidity
!
do j = 1, jmax
do i = 0, imax-1
xyr_MomFluxXCor( i,j,0 ) = xyr_MomFluxX( i,j,0 )
xyr_MomFluxYCor( i,j,0 ) = xyr_MomFluxY( i,j,0 )
xyr_HeatFluxCor( i,j,0 ) = xyr_HeatFlux( i,j,0 )
end do
end do
n = IndexH2OVap
do j = 1, jmax
do i = 0, imax-1
xyrf_QMixFluxCor( i,j,0,n ) = xyrf_QMixFlux( i,j,0,n )
end do
end do
do n = 1, IndexH2OVap-1
xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
end do
do n = IndexH2OVap+1, ncmax
xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
end do
n = IndexH2OVap
do j = 1, jmax
do i = 0, imax-1
!!$ xyr_LatentHeatFluxCor( i,j,0 ) = xy_SurfLatentHeatFlux( i,j )
xyr_LatentHeatFluxCor( i,j,0 ) = LatentHeat * xyrf_QMixFluxCor( i,j,0,n )
end do
end do
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'TauXB' , xyr_MomFluxXCor (:,:,0) )
call HistoryAutoPut( TimeN, 'TauYB' , xyr_MomFluxYCor (:,:,0) )
call HistoryAutoPut( TimeN, 'SensB' , xyr_HeatFluxCor (:,:,0) )
call HistoryAutoPut( TimeN, 'SurfH2OVapFluxB', xyrf_QMixFluxCor(:,:,0,IndexH2OVap) )
call HistoryAutoPut( TimeN, 'EvapB' , xyr_LatentHeatFluxCor(:,:,0) )
! Output of fluxes at t + \Delta t
!
! 風速, 温度, 比湿フラックス補正
! Correct fluxes of wind, temperature, specific humidity
!
do j = 1, jmax
do i = 0, imax-1
xyr_MomFluxXCor( i,j,0 ) = xyr_MomFluxX( i,j,0 ) - xy_SurfVelTransCoef( i,j ) * xyz_DUDt( i,j,1 ) * 2.0_DP * DelTime
xyr_MomFluxYCor( i,j,0 ) = xyr_MomFluxY( i,j,0 ) - xy_SurfVelTransCoef( i,j ) * xyz_DVDt( i,j,1 ) * 2.0_DP * DelTime
xyr_HeatFluxCor( i,j,0 ) = xyr_HeatFlux( i,j,0 ) - CpDry * xyr_Exner( i,j,0 ) * xy_SurfTempTransCoef( i,j ) * ( xyz_DTempDt( i,j,1 ) / xyz_Exner( i,j,1 ) - xy_DSurfTempDt( i,j ) / xyr_Exner( i,j,0 ) ) * 2.0_DP * DelTime
end do
end do
n = IndexH2OVap
do j = 1, jmax
do i = 0, imax-1
xyrf_QMixFluxCor( i,j,0,n ) = xyrf_QMixFlux( i,j,0,n ) - xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j ) * ( xyzf_DQMixDt( i,j,1,n ) - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * 2.0_DP * DelTime
end do
end do
do n = 1, IndexH2OVap-1
xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
end do
do n = IndexH2OVap+1, ncmax
xyrf_QMixFluxCor(:,:,0,n) = xyrf_QMixFlux(:,:,0,n)
end do
n = IndexH2OVap
do j = 1, jmax
do i = 0, imax-1
!!$ xyr_LatentHeatFluxCor( i,j,0 ) = xy_SurfLatentHeatFlux( i,j ) &
!!$ & - LatentHeat * xy_SurfHumidCoef( i,j ) * xy_SurfQVapTransCoef( i,j ) &
!!$ & * ( xyzf_DQMixDt( i,j,1,n ) &
!!$ & - xy_SurfDQVapSatDTemp( i,j ) * xy_DSurfTempDt( i,j ) ) * 2.0d0 * DelTime
xyr_LatentHeatFluxCor( i,j,0 ) = LatentHeat * xyrf_QMixFluxCor( i,j,0,n )
end do
end do
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'TauXA' , xyr_MomFluxXCor (:,:,0) )
call HistoryAutoPut( TimeN, 'TauYA' , xyr_MomFluxYCor (:,:,0) )
call HistoryAutoPut( TimeN, 'SensA' , xyr_HeatFluxCor (:,:,0) )
call HistoryAutoPut( TimeN, 'SurfH2OVapFluxA', xyrf_QMixFluxCor(:,:,0,IndexH2OVap) )
call HistoryAutoPut( TimeN, 'EvapA' , xyr_LatentHeatFluxCor(:,:,0) )
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'SurfH2OVapFluxU', xy_SurfH2OVapFluxA )
call HistoryAutoPut( TimeN, 'EvapU' , xy_SurfLatentHeatFluxA )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine SurfaceFluxOutput
| Subroutine : | |||
| Gamma : | real(DP), intent(in ) | ||
| xy_ZetaM(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_ZetaH(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_PsiM(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_PsiH(0:imax-1, 1:jmax) : | real(DP), intent(out)
|
Calculation of Psi_M and Psi_H
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_DP
ParamXH = ( 1.0_DP - Gamma * ZetaH )**0.25_DP
! Eq. (25)
xy_PsiM(i,j) = log( ( 1.0_DP + ParamXM )**2 * ( 1.0_DP + ParamXM**2 ) / 8.0_DP ) - 2.0_DP * atan( ParamXM ) + PI / 2.0_DP
! Eq. (26)
xy_PsiH(i,j) = log( ( 1.0_DP + ParamXH**2 )**2 / 4.0_DP )
else
! for stable condition
! Eq. (28)
xy_PsiM(i,j) = - ConstA * ZetaM - ConstB * ( ZetaM - ConstC / ConstD ) * exp( - ConstD * ZetaM ) - ConstB * ConstC / ConstD
! Eq. (32)
xy_PsiH(i,j) = - ( 1.0_DP + 2.0_DP / 3.0_DP * ConstA * ZetaH )**1.5 - ConstB * ( ZetaH - ConstC / ConstD ) * exp( - ConstD * ZetaH ) - ConstB * ConstC / ConstD + 1.0_DP
end if
end do
end do
end subroutine BH91CalcPsi
| Subroutine : | |||
| IDBulkCoefMethod : | integer , intent(in) | ||
| xy_SurfRoughLengthMom(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHeight(0:imax-1,1:jmax) : | real(DP), intent(in)
| ||
| xyz_Height(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_U(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_V(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfVirTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_VirTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfExner(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_Exner(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfVelBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfTempBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfQVapBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_BetaW(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfMOLength(0:imax-1, 1:jmax) : | real(DP), intent(out)
|
バルク係数を算出します.
Bulk coefficients are calculated.
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, xy_SurfMOLength )
!
! バルク係数を算出します.
!
! 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)
real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)
!
! Monin-Obukov length
! 作業変数
! Work variables
!
real(DP):: xy_SurfBulkRiNum (0:imax-1, 1:jmax)
! バルク $ R_i $ 数.
! Bulk $ R_i $ number
real(DP):: xy_SurfVelAbs (0:imax-1, 1:jmax)
! 風速絶対値.
! Absolute velocity
real(DP) :: xy_SurfBulkCoefMomInNeutCond (0:imax-1, 1:jmax)
real(DP) :: xy_SurfBulkCoefHeatInNeutCond (0:imax-1, 1:jmax)
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
! 実行文 ; Executable statement
!
if ( FlagConstBulkCoef ) then
! Use of constant bulk coefficient
!
xy_SurfVelBulkCoef = ConstBulkCoef
xy_SurfTempBulkCoef = ConstBulkCoef
xy_SurfQVapBulkCoef = ConstBulkCoef
xy_BetaW = 0.0_DP
xy_SurfBulkRiNum = 0.0_DP
else
select case ( IDBulkCoefMethod )
case ( IDBulkCoefMethodL82 )
! Parameterization by Louis et al. (1982)
!
! 中立バルク係数の計算
! Calculate bulk coefficient in neutral condition
!
xy_SurfBulkCoefMomInNeutCond = ( FKarm / log ( ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthMom ) / xy_SurfRoughLengthMom ) )**2
xy_SurfBulkCoefHeatInNeutCond = ( FKarm / log ( ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthMom ) / xy_SurfRoughLengthMom ) ) * ( FKarm / log ( ( xyz_Height(:,:,1) - xy_SurfHeight + xy_SurfRoughLengthHeat ) / xy_SurfRoughLengthHeat ) )
if ( FlagUseOfBulkCoefInNeutralCond ) then
! 中立条件でのバルク係数の設定
! Set bulk coefficient in neutral condition
!
xy_SurfVelBulkCoef = xy_SurfBulkCoefMomInNeutCond
xy_SurfTempBulkCoef = xy_SurfBulkCoefHeatInNeutCond
xy_SurfQVapBulkCoef = xy_SurfTempBulkCoef
xy_SurfBulkRiNum = 0.0_DP
else
! バルク $ R_i $ 数算出
! Calculate bulk $ R_i $
!
xy_SurfVelAbs = sqrt ( xy_U**2 + xy_V**2 )
xy_SurfBulkRiNum = Grav / ( xy_SurfVirTemp / xy_SurfExner ) * ( xy_VirTemp / xy_Exner - xy_SurfVirTemp / xy_SurfExner ) / max( xy_SurfVelAbs, VelMinForRi )**2 * ( xyz_Height(:,:,1) - xy_SurfHeight )
! 非中立条件でのバルク係数の計算
! Calculate bulk coefficients in non-neutral condition
!
call BulkCoefL82( xy_SurfBulkRiNum, xy_SurfRoughLengthMom , xy_SurfRoughLengthHeat, xy_SurfHeight, xyz_Height, xy_SurfBulkCoefMomInNeutCond, xy_SurfBulkCoefHeatInNeutCond, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef, xy_SurfMOLength )
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, xy_SurfMOLength )
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)
| ||
| xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHeight(0:imax-1,1:jmax) : | real(DP), intent(in)
| ||
| xyz_Height(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_U(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_V(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfVirTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_VirTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfExner(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_Exner(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfVelBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfTempBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfQVapBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_BetaW(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfBulkRiNum(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfMOLength(0:imax-1, 1:jmax) : | real(DP), intent(out) |
バルク係数を算出します.
Bulk coefficients are calculated.
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, xy_SurfMOLength )
!
! バルク係数を算出します.
!
! 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
real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)
! 作業変数
! Work variables
!
real(DP):: xy_SurfVelAbs (0:imax-1, 1:jmax)
! 風速絶対値.
! Absolute velocity
real(DP) :: xy_MOLength (0:imax-1, 1:jmax)
real(DP) :: xy_MOLengthSave(0:imax-1, 1:jmax)
real(DP) :: xy_ZetaM(0:imax-1, 1:jmax)
real(DP) :: xy_ZetaH(0:imax-1, 1:jmax)
real(DP) :: xy_PsiM1(0:imax-1, 1:jmax)
real(DP) :: xy_PsiH1(0:imax-1, 1:jmax)
real(DP) :: xy_PsiM0(0:imax-1, 1:jmax)
real(DP) :: xy_PsiH0(0:imax-1, 1:jmax)
real(DP) :: xy_FricVelByU1(0:imax-1, 1:jmax)
real(DP) :: SurfBulkRiNum
logical :: FlagConverge
logical :: xy_FlagConverge(0:imax-1, 1:jmax)
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.0e10_DP
xy_BetaW = 0.0_DP
xy_FlagConverge = .false.
loop_iteration : do iLoop = 1, nLoop
xy_MOLengthSave = xy_MOLength
! Calculation of Psi_{M,H}
if ( FlagUseOfBulkCoefInNeutralCond ) then
xy_PsiM0 = 0.0_DP
xy_PsiH0 = 0.0_DP
xy_PsiM1 = 0.0_DP
xy_PsiH1 = 0.0_DP
else
xy_ZetaM = xy_SurfRoughLengthMom / xy_MOLength
xy_ZetaH = xy_SurfRoughLengthHeat / xy_MOLength
call BH91CalcPsi( Gamma, 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.0e-10_DP
xy_MOLength(i,j) = ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) ) / ( FKarm * SurfBulkRiNum ) * xy_SurfVelBulkCoef(i,j)**1.5_DP / xy_SurfTempBulkCoef(i,j)
end if
end do
end do
! TO BE DELETED
!!$ ! Check of convergence
!!$ FlagConverge = .true.
!!$ loop_check : do j = 1, jmax
!!$ do i = 0, imax-1
!!$ MOLengthErr = &
!!$ & abs( xy_MOLength(i,j) - xy_MOLengthSave(i,j) ) &
!!$ & / max( abs( xy_MOLength(i,j) ), 1.0d-10 )
!!$ if ( MOLengthErr > MOLengthErrCriterion ) then
!!$ FlagConverge = .false.
!!$ exit loop_check
!!$ end if
!!$ end do
!!$ end do loop_check
!!$ a_FlagReCalcLocal = ( .not. FlagConverge )
!!$ call MPIWrapperChkTrue( &
!!$ & 1, a_FlagReCalcLocal, & ! (in)
!!$ & a_FlagReCalcGlobal & ! (out)
!!$ & )
!!$ if ( .not. a_FlagReCalcGlobal(1) ) exit loop_iteration
! Check of convergence
do j = 1, jmax
do i = 0, imax-1
if ( .not. xy_FlagConverge(i,j) ) then
MOLengthErr = abs( xy_MOLength(i,j) - xy_MOLengthSave(i,j) ) / max( abs( xy_MOLength(i,j) ), 1.0e-10_DP )
if ( MOLengthErr <= MOLengthErrCriterion ) then
xy_FlagConverge(i,j) = .true.
end if
end if
end do
end do
FlagConverge = .true.
loop_check : do j = 1, jmax
do i = 0, imax-1
if ( .not. xy_FlagConverge(i,j) ) then
FlagConverge = .false.
exit loop_check
end if
end do
end do loop_check
if ( FlagConverge ) exit loop_iteration
end do loop_iteration
if ( iLoop > nLoop ) then
call MessageNotify( 'E', module_name, 'Monin-Obukhov length is not convergent.' )
end if
xy_SurfMOLength = xy_MOLength
call HistoryAutoPut( TimeN, 'MOLength' , xy_MOLength )
call HistoryAutoPut( TimeN, 'MOLengthInv', 1.0_DP/xy_MOLength )
call HistoryAutoPut( TimeN, 'BetaW' , xy_BetaW )
call HistoryAutoPut( TimeN, 'BLHeight' , xy_BLHeight )
end subroutine BulkCoefBH91B94
| Subroutine : | |||
| xy_SurfBulkRiNum(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfRoughLengthMom(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfRoughLengthHeat(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHeight(0:imax-1,1:jmax) : | real(DP), intent(in)
| ||
| xyz_Height(0:imax-1, 1:jmax, 1:kmax) : | 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(in) | ||
| xy_SurfVelBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfTempBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfQVapBulkCoef(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfMOLength(0:imax-1, 1:jmax) : | real(DP), intent(out) |
バルク係数を算出します.
Bulk coefficients are calculated.
subroutine BulkCoefL82( xy_SurfBulkRiNum, xy_SurfRoughLengthMom , xy_SurfRoughLengthHeat, xy_SurfHeight, xyz_Height, xy_SurfBulkCoefMomInNeutCond, xy_SurfBulkCoefHeatInNeutCond, xy_SurfVelBulkCoef, xy_SurfTempBulkCoef, xy_SurfQVapBulkCoef, xy_SurfMOLength )
!
! バルク係数を算出します.
!
! Bulk coefficients are calculated.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm ! $ k $ .
! カルマン定数.
! Karman constant
! 宣言文 ; Declaration statements
!
real(DP), intent(in):: xy_SurfBulkRiNum (0:imax-1, 1:jmax)
! バルク $ R_i $ 数.
! Bulk $ R_i $ number
real(DP), intent(in):: xy_SurfRoughLengthMom (0:imax-1, 1:jmax)
! 地表粗度長
! Surface rough length for momentum
real(DP), intent(in):: xy_SurfRoughLengthHeat(0:imax-1, 1:jmax)
! 地表粗度長
! Surface rough length for heat
real(DP), intent(in):: xy_SurfHeight(0:imax-1,1:jmax)
! $ z_s $ . 地表面高度.
! Surface height.
real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
! 高度.
! Height
real(DP), intent(in):: xy_SurfBulkCoefMomInNeutCond (0:imax-1, 1:jmax)
!
!
real(DP), intent(in):: xy_SurfBulkCoefHeatInNeutCond(0:imax-1, 1:jmax)
!
!
real(DP), intent(out):: xy_SurfVelBulkCoef (0:imax-1, 1:jmax)
! バルク係数:運動量.
! Bulk coefficient: temperature
real(DP), intent(out):: xy_SurfTempBulkCoef (0:imax-1, 1:jmax)
! バルク係数:温度.
! Bulk coefficient: temperature
real(DP), intent(out):: xy_SurfQVapBulkCoef (0:imax-1, 1:jmax)
! バルク係数:比湿.
! Bulk coefficient: specific humidity
real(DP), intent(out):: xy_SurfMOLength(0:imax-1, 1:jmax)
! 作業変数
! Work variables
!
real(DP) :: SurfBulkRiNum
real(DP) :: xy_MOLength(0:imax-1, 1:jmax)
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
! 実行文 ; Executable statement
!
! 非中立条件でのバルク係数の計算
! Calculate bulk coefficients in non-neutral condition
!
! Parameterization by Louis et al. (1982)
!
do j = 1, jmax
do i = 0, imax-1
if ( xy_SurfBulkRiNum(i,j) > 0.0_DP ) then
xy_SurfVelBulkCoef(i,j) = xy_SurfBulkCoefMomInNeutCond(i,j) / ( 1.0_DP + 10.0_DP * xy_SurfBulkRiNum(i,j) / sqrt( 1.0_DP + 5.0_DP * xy_SurfBulkRiNum(i,j) ) )
xy_SurfTempBulkCoef(i,j) = xy_SurfBulkCoefHeatInNeutCond(i,j) / ( 1.0_DP + 15.0_DP * xy_SurfBulkRiNum(i,j) * sqrt( 1.0_DP + 5.0_DP * xy_SurfBulkRiNum(i,j) ) )
xy_SurfQVapBulkCoef(i,j) = xy_SurfTempBulkCoef(i,j)
else
xy_SurfVelBulkCoef(i,j) = xy_SurfBulkCoefMomInNeutCond(i,j) * ( 1.0_DP - 10.0_DP * xy_SurfBulkRiNum(i,j) / ( 1.0_DP + 75.0_DP * xy_SurfBulkCoefMomInNeutCond(i,j) * sqrt( - ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) + xy_SurfRoughLengthMom(i,j) ) / xy_SurfRoughLengthMom(i,j) * xy_SurfBulkRiNum(i,j) ) ) )
xy_SurfTempBulkCoef(i,j) = xy_SurfBulkCoefHeatInNeutCond(i,j) * ( 1.0_DP - 15.0_DP * xy_SurfBulkRiNum(i,j) / ( 1.0_DP + 75.0_DP * xy_SurfBulkCoefHeatInNeutCond(i,j) * sqrt( - ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) + xy_SurfRoughLengthHeat(i,j) ) / xy_SurfRoughLengthHeat(i,j) * xy_SurfBulkRiNum(i,j) ) ) )
xy_SurfQVapBulkCoef(i,j) = xy_SurfTempBulkCoef(i,j)
end if
! Calculation of Monin-Obukhov length
SurfBulkRiNum = xy_SurfBulkRiNum(i,j)
if ( SurfBulkRiNum == 0.0_DP ) SurfBulkRiNum = 1.0e-10_DP
xy_MOLength(i,j) = ( xyz_Height(i,j,1) - xy_SurfHeight(i,j) ) / ( FKarm * SurfBulkRiNum ) * xy_SurfVelBulkCoef(i,j)**1.5_DP / xy_SurfTempBulkCoef(i,j)
end do
end do
xy_SurfMOLength = xy_MOLength
end subroutine BulkCoefL82
| Variable : | |||
| FlagConstBulkCoef : | logical, save
|
| Variable : | |||
| FlagUseOfBulkCoefInNeutralCond : | logical, save
|
| Variable : | |||
| FricTimeConstAtLB : | real(DP), save
|
| Variable : | |||
| HeatFluxAtLB : | real(DP), save
|
| Variable : | |||
| LowLatFricAtLB : | real(DP), save
|
| Variable : | |||
| MassFluxAtLB : | real(DP), save
|
| Variable : | |||
| QVapBulkCoefMax : | real(DP), save
|
| Variable : | |||
| QVapBulkCoefMin : | real(DP), save
|
| Variable : | |||
| TempBulkCoefMax : | real(DP), save
|
| Variable : | |||
| TempBulkCoefMin : | real(DP), save
|
| Variable : | |||
| VelBulkCoefMax : | real(DP), save
|
| Variable : | |||
| VelBulkCoefMin : | real(DP), save
|
| Variable : | |||
| VelMaxForQVap : | real(DP), save
|
| Variable : | |||
| VelMaxForTemp : | real(DP), save
|
| Variable : | |||
| VelMaxForVel : | real(DP), save
|
| Variable : | |||
| VelMinForQVap : | real(DP), save
|
| Variable : | |||
| VelMinForRi : | real(DP), save
|
| Variable : | |||
| VelMinForTemp : | real(DP), save
|
| Variable : | |||
| VelMinForVel : | real(DP), save
|
| Constant : | |||
| module_name = ‘surface_flux_bulk‘ : | character(*), parameter
|
| Variable : | |||
| surface_flux_bulk_inited = .false. : | logical, save
|
| Constant : | |||
| version = ’$Name: $’ // ’$Id: surface_flux_bulk.f90,v 1.27 2015/02/06 11:25:14 yot Exp $’ : | character(*), parameter
|