| Class | vdiffusion_my |
| In: |
vdiffusion/vdiffusion_my.f90
|
Note that Japanese and English are described in parallel.
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated.
| VDiffusion : | 鉛直拡散フラックスの計算 |
| VDiffusionOutPut : | フラックスの出力 |
| ———— : | ———— |
| VDiffusion : | Calculate vertical diffusion fluxes |
| VDiffusionOutPut : | Output fluxes |
| Subroutine : | |||
| 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)
| ||
| xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | 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)
| ||
| 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)
| ||
| xyr_Height(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)
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated by use of MY2 model.
subroutine VDiffusion( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated by use of MY2 model.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm, Grav
! $ g $ [m s-2].
! 重力加速度.
! Gravitational acceleration
! 時刻管理
! Time control
!
use timeset, only: TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! $ u $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! $ v $ . 南北風速. Northward wind
real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q $ . 質量混合比. Mass mixing ratio
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T} $ . 温度 (半整数レベル).
! Temperature (half level)
real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
! $ T_v $ . 仮温度. Virtual temperature
real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air 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):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
! 高度 (半整数レベル).
! Height (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(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 質量フラックス.
! Mass flux of compositions
real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
! 作業変数
! Work variables
!
real(DP) :: xyr_DVelDz (0:imax-1, 1:jmax, 0:kmax)
! $ \DD{|\Dvect{v}|}{z} $
real(DP) :: xyr_BulkRiNum (0:imax-1, 1:jmax, 0:kmax)
! バルク $ R_i $ 数.
! Bulk $ R_i $
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! 拡散係数の計算
! Calculation of diffusion coefficient
!
if ( FlagConstDiffCoef ) then
xyr_VelDiffCoef (:,:,0 ) = 0.0_DP
xyr_VelDiffCoef (:,:,1:kmax-1) = ConstDiffCoefM
xyr_VelDiffCoef (:,:,kmax ) = 0.0_DP
xyr_TempDiffCoef(:,:,0 ) = 0.0_DP
xyr_TempDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH
xyr_TempDiffCoef(:,:,kmax ) = 0.0_DP
xyr_QMixDiffCoef(:,:,0 ) = 0.0_DP
xyr_QMixDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH
xyr_QMixDiffCoef(:,:,kmax ) = 0.0_DP
else
! バルク $ R_i $ 数算出
! Calculate bulk $ R_i $
!
xyr_DVelDz(:,:,0) = 0.0_DP
xyr_DVelDz(:,:,kmax) = 0.0_DP
xyr_BulkRiNum(:,:,0) = 0.0_DP
xyr_BulkRiNum(:,:,kmax) = 0.0_DP
do k = 1, kmax-1
xyr_DVelDz(:,:,k) = sqrt( max( SquareVelMin , ( xyz_U(:,:,k+1) - xyz_U(:,:,k) )**2 + ( xyz_V(:,:,k+1) - xyz_V(:,:,k) )**2 ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
xyr_BulkRiNum(:,:,k) = Grav / ( xyr_VirTemp(:,:,k) / xyr_Exner(:,:,k) ) * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) / xyr_DVelDz(:,:,k)**2
xyr_BulkRiNum(:,:,k) = max( xyr_BulkRiNum(:,:,k) , BulkRiNumMin )
end do
! 拡散係数の計算
! Calculate diffusion coefficients
!
call VDiffCoefficient( xy_SurfHeight, xyr_Height, xyr_DVelDz, xyr_BulkRiNum, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )
end if
! 浅い積雲対流
! Shallow cumulus convection
!
! (AGCM5 から導入予定)
! 拡散係数の出力
! Output diffusion coefficients
!
! (上記の「浅い積雲対流」導入後に作成)
! 拡散係数出力
! Diffusion coeffficients output
!
call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef )
call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef )
call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef )
! 輸送係数とフラックスの計算
! Calculate transfer coefficient and flux
!
call VDiffusionCalcFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine VDiffusion
| Subroutine : | |||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ), optional
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ), optional
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ), optional
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(in ), optional
| ||
| xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out), optional
| ||
| xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out), optional
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out), optional
| ||
| xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out), optional
|
時間変化率の計算を行います.
Calculate tendencies.
subroutine VDiffusionExpTendency( xyr_Press, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt )
!
! 時間変化率の計算を行います.
!
! Calculate tendencies.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in ), optional :: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(in ), optional :: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(in ), optional :: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in ), optional :: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 比湿フラックス.
! Specific humidity flux
real(DP), intent(out), optional :: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{u}{t} $ . 東西風速変化.
! Eastward wind tendency
real(DP), intent(out), optional :: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{v}{t} $ . 南北風速変化.
! Northward wind tendency
real(DP), intent(out), optional :: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(out), optional :: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 質量混合比変化.
! Mass mixing ratio tendency
! 作業変数
! Work variables
!
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! Check arguments
!
if ( present( xyz_DUDt ) ) then
if ( .not. present( xyr_MomFluxX ) ) then
call MessageNotify( 'E', module_name, 'xyr_MomFluxX has to be present.' )
end if
end if
if ( present( xyz_DVDt ) ) then
if ( .not. present( xyr_MomFluxY ) ) then
call MessageNotify( 'E', module_name, 'xyr_MomFluxY has to be present.' )
end if
end if
if ( present( xyz_DTempDt ) ) then
if ( .not. present( xyr_HeatFlux ) ) then
call MessageNotify( 'E', module_name, 'xyr_HeatFlux has to be present.' )
end if
end if
if ( present( xyzf_DQMixDt ) ) then
if ( .not. present( xyrf_QMixFlux ) ) then
call MessageNotify( 'E', module_name, 'xyrf_QMixFlux has to be present.' )
end if
end if
if ( present( xyz_DUDt ) ) then
do k = 1, kmax
xyz_DUDt(:,:,k) = + Grav * ( xyr_MomFluxX(:,:,k) - xyr_MomFluxX(:,:,k-1) ) / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) )
end do
end if
if ( present( xyz_DVDt ) ) then
do k = 1, kmax
xyz_DVDt(:,:,k) = + Grav * ( xyr_MomFluxY(:,:,k) - xyr_MomFluxY(:,:,k-1) ) / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) )
end do
end if
if ( present( xyz_DTempDt ) ) then
do k = 1, kmax
xyz_DTempDt(:,:,k) = + Grav / CpDry * ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) ) / ( xyr_Press (:,:,k) - xyr_Press (:,:,k-1) )
end do
end if
if ( present( xyzf_DQMixDt ) ) then
do n = 1, ncmax
do k = 1, kmax
xyzf_DQMixDt(:,:,k,n) = + Grav * ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) ) / ( xyr_Press (:,:,k) - xyr_Press (:,:,k-1) )
end do
end do
end if
end subroutine VDiffusionExpTendency
| Subroutine : |
vdiffusion_my モジュールの初期化を行います. NAMELIST#vdiffusion_my_nml の読み込みはこの手続きで行われます.
"vdiffusion_my" module is initialized. "NAMELIST#vdiffusion_my_nml" is loaded in this procedure.
This procedure input/output NAMELIST#vdiffusion_my_nml .
subroutine VDiffusionInit
!
! vdiffusion_my モジュールの初期化を行います.
! NAMELIST#vdiffusion_my_nml の読み込みはこの手続きで行われます.
!
! "vdiffusion_my" module is initialized.
! "NAMELIST#vdiffusion_my_nml" is loaded in this procedure.
!
! モジュール引用 ; USE statements
!
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output
! 文字列操作
! Character handling
!
use dc_string, only: StoA
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoAddVariable
! 座標データ設定
! Axes data settings
!
use axesset, only: AxnameX, AxnameY, AxnameZ, AxnameR, AxnameT
! 宣言文 ; Declaration statements
!
implicit none
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
namelist /vdiffusion_my_nml/ FlagConstDiffCoef, ConstDiffCoefM, ConstDiffCoefH, SquareVelMin, BulkRiNumMin, MixLengthMax, ShMin, SmMin, VelDiffCoefMin, TempDiffCoefMin, VelDiffCoefMax, TempDiffCoefMax, MYConstA1, MYConstB1, MYConstA2, MYConstB2, MYConstC1
!
! デフォルト値については初期化手続 "vdiffusion_my#VDiffInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "vdiffusion_my#VDiffInit" for the default values.
!
! 実行文 ; Executable statement
!
if ( vdiffusion_my_inited ) return
! デフォルト値の設定
! Default values settings
!
FlagConstDiffCoef = .false.
ConstDiffCoefM = 0.0_DP
ConstDiffCoefH = 0.0_DP
SquareVelMin = 0.1_DP
BulkRiNumMin = - 100.0_DP
MixLengthMax = 300.0_DP
ShMin = 0.0_DP
SmMin = 0.0_DP
VelDiffCoefMin = 0.1_DP
TempDiffCoefMin = 0.1_DP
VelDiffCoefMax = 10000.0_DP
TempDiffCoefMax = 10000.0_DP
! Parameters proposed by Mellor and Yamada (1982).
!
MYConstA1 = 0.92_DP
MYConstB1 = 16.6_DP
MYConstA2 = 0.74_DP
MYConstB2 = 10.1_DP
MYConstC1 = 0.08_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 = vdiffusion_my_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
if ( iostat_nml == 0 ) write( STDOUT, nml = vdiffusion_my_nml )
end if
! ヒストリデータ出力のためのへの変数登録
! Register of variables for history data output
!
call HistoryAutoAddVariable( 'VelDiffCoef', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'diffusion coef. momentum', 'm2 s-1' )
call HistoryAutoAddVariable( 'TempDiffCoef', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'diffusion coef. heat ', 'm2 s-1' )
call HistoryAutoAddVariable( 'QVapDiffCoef', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'diffusion coef. moisture', 'm2 s-1' )
call HistoryAutoAddVariable( 'MomFluxX', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'eastward momentum flux', 'N m-2' )
call HistoryAutoAddVariable( 'MomFluxY', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'northward momentum flux', 'N m-2' )
call HistoryAutoAddVariable( 'HeatFlux', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'heat flux', 'W m-2' )
call HistoryAutoAddVariable( 'QVapFlux', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'moisture flux', 'W m-2' )
call HistoryAutoAddVariable( 'DUDtVDiff', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'tendency of zonal wind by vertical diffusion', 'm s-2' )
call HistoryAutoAddVariable( 'DVDtVDiff', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'tendency of meridional wind by vertical diffusion', 'm s-2' )
call HistoryAutoAddVariable( 'DTempDtVDiff', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'tendency of temperature by vertical diffusion', 'K s-1' )
call HistoryAutoAddVariable( 'DQVapDtVDiff', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'tendency of specific humidity by vertical diffusion', 's-1' )
call HistoryAutoAddVariable( 'TurKinEne', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'turbulent kinetic energy', 'm2 s-2' )
call HistoryAutoAddVariable( 'TKEPShear', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'turbulent kinetic energy production rate by shear', 'm2 s-3' )
call HistoryAutoAddVariable( 'TKEPBuoy', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'turbulent kinetic energy production rate by buoyancy', 'm2 s-3' )
call HistoryAutoAddVariable( 'TKEDiss', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'turbulent kinetic energy dissipation rate', 'm2 s-3' )
call HistoryAutoAddVariable( 'MixLength', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'mixing length', 'm' )
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, 'For vertical diffusion flux:' )
call MessageNotify( 'M', module_name, ' FlagConstDiffCoef = %b', l = (/ FlagConstDiffCoef /) )
call MessageNotify( 'M', module_name, ' ConstDiffCoefM = %f', d = (/ ConstDiffCoefM /) )
call MessageNotify( 'M', module_name, ' ConstDiffCoefH = %f', d = (/ ConstDiffCoefH /) )
call MessageNotify( 'M', module_name, ' SquareVelMin = %f', d = (/ SquareVelMin /) )
call MessageNotify( 'M', module_name, ' BulkRiNumMin = %f', d = (/ BulkRiNumMin /) )
call MessageNotify( 'M', module_name, 'For diffusion coefficients:' )
call MessageNotify( 'M', module_name, ' MixLengthMax = %f', d = (/ MixLengthMax /) )
call MessageNotify( 'M', module_name, ' ShMin = %f', d = (/ ShMin /) )
call MessageNotify( 'M', module_name, ' SmMin = %f', d = (/ SmMin /) )
call MessageNotify( 'M', module_name, ' VelDiffCoefMin = %f', d = (/ VelDiffCoefMin /) )
call MessageNotify( 'M', module_name, ' TempDiffCoefMin = %f', d = (/ TempDiffCoefMin /) )
call MessageNotify( 'M', module_name, ' VelDiffCoefMax = %f', d = (/ VelDiffCoefMax /) )
call MessageNotify( 'M', module_name, ' TempDiffCoefMax = %f', d = (/ TempDiffCoefMax /) )
call MessageNotify( 'M', module_name, ' MYConstA1 = %f', d = (/ MYConstA1 /) )
call MessageNotify( 'M', module_name, ' MYConstB1 = %f', d = (/ MYConstB1 /) )
call MessageNotify( 'M', module_name, ' MYConstA2 = %f', d = (/ MYConstA2 /) )
call MessageNotify( 'M', module_name, ' MYConstB2 = %f', d = (/ MYConstB2 /) )
call MessageNotify( 'M', module_name, ' MYConstC1 = %f', d = (/ MYConstC1 /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
vdiffusion_my_inited = .true.
end subroutine VDiffusionInit
| Subroutine : | |||
| 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)
| ||
| xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | 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)
| ||
| 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)
| ||
| xyr_Height(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)
| ||
| xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_SurfMomFluxX(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfMomFluxY(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyz_DTurKinEneDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated by use of MY2.5 model.
subroutine VDiffusionMY25( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyz_TurKinEne, xy_SurfMomFluxX, xy_SurfMomFluxY, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyz_DTurKinEneDt )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated by use of MY2.5 model.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm, Grav, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! MPI 関連ルーチン
! MPI related routines
!
use mpi_wrapper, only : MPIWrapperChkTrue
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! $ u $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! $ v $ . 南北風速. Northward wind
real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q $ . 質量混合比. Mass mixing ratio
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T} $ . 温度 (半整数レベル).
! Temperature (half level)
real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
! $ T_v $ . 仮温度. Virtual temperature
real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air 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):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
! 高度 (半整数レベル).
! Height (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):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax)
!
! Turbulent kinetic energy (m2 s-2)
real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax)
!
! Eastward momentum flux at surface
real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax)
!
! Northward momentum flux at surface
real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 質量フラックス.
! Mass flux of compositions
real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Diffusion coefficient: temperature
real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax)
!
! Tendency of turbulent kinetic energy
! 作業変数
! Work variables
!
real(DP) :: xyz_MixLength(0:imax-1, 1:jmax, 1:kmax)
! 混合距離.
! Mixing length
real(DP) :: xyz_DVelDzSq(0:imax-1, 1:jmax, 1:kmax)
!
! Vertical shear squared (s-2)
real(DP) :: xyz_StatStab(0:imax-1, 1:jmax, 1:kmax)
!
! Static stability (s-2)
real(DP) :: GhMin
!
! Minimum of G_h
real(DP) :: GhMax
!
! Maximum of G_h
real(DP) :: xyz_Gm(0:imax-1, 1:jmax, 1:kmax)
!
! G_m
real(DP) :: xyz_Gh(0:imax-1, 1:jmax, 1:kmax)
!
! G_h
!!$ real(DP) :: xyz_SmHat(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! \hat{S}_M
!!$ real(DP) :: xyz_ShHat(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! \hat{S}_h
!!$ real(DP) :: xyz_DSmHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! derivative of \hat{S}_M
!!$ real(DP) :: xyz_DShHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! derivative of \hat{S}_h
real(DP) :: xyz_Sm(0:imax-1, 1:jmax, 1:kmax)
!
! S_M
real(DP) :: xyz_Sh(0:imax-1, 1:jmax, 1:kmax)
!
! S_h
real(DP), parameter :: Stke = 0.2_DP
!
! S_{TKE} = 0.2
real(DP) :: xyz_VelDiffCoef (0:imax-1, 1:jmax, 1:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP) :: xyz_TempDiffCoef (0:imax-1, 1:jmax, 1:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_TurKinEneDiffCoef (0:imax-1, 1:jmax, 0:kmax)
!
! Diffusion coefficient: turbulent kinetic energy
real(DP) :: xyz_TurKinEneDiffCoef (0:imax-1, 1:jmax, 1:kmax)
!
! Diffusion coefficient: turbulent kinetic energy
real(DP) :: xyr_TurKinEneTransCoef(0:imax-1, 1:jmax, 0:kmax)
!
! Transfer coefficient: turbulent kinetic energy
real(DP) :: xyr_TurKinEneFlux(0:imax-1, 1:jmax, 0:kmax)
!
! Turbulent energy flux
real(DP) :: xyz_CShe1(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CShe2(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CBuo1(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CBuo2(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CDis1(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CDis2(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_TurKinEneProShear(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_TurKinEneProBuoya(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_FricVelSq (0:imax-1, 1:jmax)
real(DP) :: xy_TurKinEneAtLB(0:imax-1, 1:jmax)
real(DP) :: xyza_TurKinEneMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
!
! Implicit matrix for turbulent kinetic energy
real(DP) :: xyz_TurKinEneVec(0:imax-1, 1:jmax, 1:kmax)
!
! Implicit vector for turbulent kinetic energy
real(DP) :: xyza_TurKinEneLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP) :: xyz_DelTurKinEneLUVec(0:imax-1, 1:jmax, 1:kmax)
!
! Tendency of turbulent kinetic energy
real(DP) :: xyz_TurKinEneDiss(0:imax-1, 1:jmax, 1:kmax)
!
! Dissipation rate of turbulent kinetic energy (m2 s-3)
real(DP) :: xyz_TurKinEneNonZero(0:imax-1, 1:jmax, 1:kmax)
!
! Turbulent kinetic energy with offset (m2 s-2)
real(DP), parameter :: TurKinEneOffset = ( 1.0e-3_DP )**2 / 2.0_DP
logical :: FlagReCalc
!
! Flag for recalculation
logical :: a_FlagReCalcLocal (1)
logical :: a_FlagReCalcGlobal(1)
integer :: iloop
integer :: nloop
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_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 turbulent kinetic energy with offset
!
xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset
!
! Calculation of vertical shear squared
do k = 1, kmax
if ( k == 1 ) then
xyz_DVelDzSq(:,:,k) = ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k ) )**2 + ( xyz_V(:,:,k+1) - xyz_V(:,:,k ) )**2 ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) )**2
!!$ xyz_DVelDzSq(:,:,k) = &
!!$ & ( ( xyz_U(:,:,k+1) - 0.0_DP )**2 &
!!$ & + ( xyz_V(:,:,k+1) - 0.0_DP )**2 ) &
!!$ & / ( xyz_Height(:,:,k+1) - xy_SurfHeight )**2
else if ( k == kmax ) then
xyz_DVelDzSq(:,:,k) = ( ( xyz_U(:,:,k ) - xyz_U(:,:,k-1) )**2 + ( xyz_V(:,:,k ) - xyz_V(:,:,k-1) )**2 ) / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) )**2
else
xyz_DVelDzSq(:,:,k) = ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k-1) )**2 + ( xyz_V(:,:,k+1) - xyz_V(:,:,k-1) )**2 ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )**2
end if
end do
! Calculation of static stability
do k = 1, kmax
if ( k == 1 ) then
xyz_StatStab(:,:,k) = Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) )
else if ( k == kmax ) then
xyz_StatStab(:,:,k) = Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) * ( xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) )
else
xyz_StatStab(:,:,k) = Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )
end if
end do
! 混合距離の算出
! Calculate mixing length
!
do k = 1, kmax
xyz_MixLength(:,:,k) = FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / (1.0_DP + FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / MixLengthMax )
end do
! Limit mixing length (Galperin et al., 1988) and avoid zero
xyz_MixLength = min( xyz_MixLength, 0.53_DP * sqrt( 2.0_DP * xyz_TurKinEneNonZero / max( xyz_StatStab, 1.0e-10_DP ) ) ) + 1.0e-10_DP
xyz_Gh = - xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_StatStab
! Actually, xyz_Gm is not used below.
xyz_Gm = xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_DVelDzSq
! Limit Gh (Galperin et al., 1988)
GhMin = - 0.53_DP**2
GhMax = 1.0_DP / ( MYConstA2 * ( 12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) )
xyz_Gh = max( GhMin, min( xyz_Gh, GhMax ) )
xyz_Sh = MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) / ( 1.0_DP - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_Gh )
xyz_Sm = ( MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 - 6.0_DP * MYConstA1 / MYConstB1 ) + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) * xyz_Gh * xyz_Sh ) / ( 1.0_DP - 9.0_DP * MYConstA1 * MYConstA2 * xyz_Gh )
!!$ xyz_DShHatDTKE = &
!!$ & - 2.0_DP * MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) &
!!$ & / ( 2.0_DP * xyz_TurKinEneNonZero &
!!$ & - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_GhPrime )**2
!!$ xyz_DSmHatDTKE = &
!!$ & ( &
!!$ & 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) &
!!$ & * xyz_GhPrime * xyz_DShHatDTKE &
!!$ & * ( 2.0_DP * xyz_TurKinEneNonZero &
!!$ & - 9.0_DP * MYConstA1 * MYConstA2 * xyz_GhPrime ) &
!!$ & - 2.0_DP &
!!$ & * ( &
!!$ & MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 &
!!$ & - 6.0_DP * MYConstA1 / MYConstB1 ) &
!!$ & + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) &
!!$ & * xyz_GhPrime * xyz_ShHat &
!!$ & ) &
!!$ & ) &
!!$ & / ( 2.0_DP * xyz_TurKinEneNonZero &
!!$ & - 9.0_DP * MYConstA1 * MYConstA2 * xyz_GhPrime )**2
! 拡散係数の計算
! Calculation of diffusion coefficient
!
xyz_VelDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm
xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh
!
do k = 0, kmax
if ( ( k == 0 ) .or. ( k == kmax ) ) then
xyr_VelDiffCoef (:,:,k) = 0.0_DP
xyr_TempDiffCoef(:,:,k) = 0.0_DP
else
xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP
xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP
end if
end do
!
do k = 1, kmax-1
do j = 1, jmax
do i = 0, imax-1
xyr_VelDiffCoef(i,j,k) = max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), VelDiffCoefMin )
xyr_TempDiffCoef(i,j,k) = max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), TempDiffCoefMin )
end do
end do
end do
!
xyr_QMixDiffCoef = xyr_TempDiffCoef
! 輸送係数とフラックスの計算
! Calculate transfer coefficient and flux
!
call VDiffusionCalcFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux )
! Calculate tendency of turbulent kinetic energy
! Set diffusion coefficient for turbulent kinetic energy
xyz_TurKinEneDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * Stke
!
do k = 0, kmax
if ( k == 0 ) then
xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1)
else if ( k == kmax ) then
xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,kmax)
else
xyr_TurKinEneDiffCoef(:,:,k) = ( xyz_TurKinEneDiffCoef(:,:,k) + xyz_TurKinEneDiffCoef(:,:,k+1) ) / 2.0_DP
end if
end do
! Calculate turbulent kinetic energy at lower boundary
!
xy_FricVelSq = sqrt( xy_SurfMomFluxX**2 + xy_SurfMomFluxY**2 ) / ( xyr_Press(:,:,0) / ( GasRDry * xyr_VirTemp(:,:,0) ) )
xy_TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * xy_FricVelSq
xy_TurKinEneAtLB = xy_TurKinEneAtLB + TurKinEneOffset
! Calculate transfer coefficient and flux of turbulent kinetic energy
!
! When transfer coefficient at lower boundary is calculated,
! diffusion coefficient at mid-point of 1st layer is used.
! In addition, transfer coefficient at upper boundary is assumed
! to be zero.
k = 0
xyr_TurKinEneTransCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,1) - xy_SurfHeight )
do k = 1, kmax-1
xyr_TurKinEneTransCoef(:,:,k) = xyr_TurKinEneDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
end do
k = kmax
xyr_TurKinEneTransCoef(:,:,k) = 0.0_DP
!
do k = 1, kmax-1
xyr_TurKinEneFlux(:,:,k) = - xyr_TurKinEneTransCoef(:,:,k) * ( xyz_TurKinEneNonZero(:,:,k+1) - xyz_TurKinEneNonZero(:,:,k) )
end do
k = 0
xyr_TurKinEneFlux(:,:,k) = - xyr_TurKinEneTransCoef(:,:,k) * ( xyz_TurKinEneNonZero(:,:,k+1) - xy_TurKinEneAtLB )
k = kmax
xyr_TurKinEneFlux(:,:,k) = 0.0_DP
!!$ xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_SmHat &
!!$ & * xyz_DVelDzSq &
!!$! & * xyz_TurKinEne**1.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$ xyz_CShe2 = 1.5_DP * 2.0d0**1.5 * xyz_MixLength * xyz_SmHat &
!!$ & * xyz_DVelDzSq &
!!$! & * xyz_TurKinEne**0.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$ xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_ShHat &
!!$ & * xyz_StatStab &
!!$! & * xyz_TurKinEne**1.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$ xyz_CBuo2 = - 1.5_DP * 2.0d0**1.5 * xyz_MixLength * xyz_ShHat &
!!$ & * xyz_StatStab &
!!$! & * xyz_TurKinEne**0.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$ xyz_CDis1 = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$! & * xyz_TurKinEne**1.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$ xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$! & * xyz_TurKinEne**0.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$! xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_SmHat &
!!$! & * xyz_DVelDzSq &
!!$!! & * xyz_TurKinEne**1.5
!!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$ xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_Sm &
!!$ & * xyz_DVelDzSq &
!!$! & * xyz_TurKinEne**1.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$ xyz_CShe2 = 0.0_DP
!!$! xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_ShHat &
!!$! & * xyz_StatStab &
!!$!! & * xyz_TurKinEne**1.5
!!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$ xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_Sh &
!!$ & * xyz_StatStab &
!!$! & * xyz_TurKinEne**1.5
!!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$ xyz_CBuo2 = 0.0_DP
!!$ xyz_CDis1 = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$ & * xyz_TurKinEne**1.5
!!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$ xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$ & * xyz_TurKinEne**0.5
!!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$ xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_DVelDzSq &
!!$ & * xyz_SmHat &
!!$ & * xyz_TurKinEneNonZero**1.5
!!$ xyz_CShe2 = 2.0d0**1.5 * xyz_MixLength * xyz_DVelDzSq &
!!$ & * ( xyz_DSmHatDTKE &
!!$ & * xyz_TurKinEneNonZero**1.5 &
!!$ & + 1.5_DP * xyz_SmHat &
!!$ & * xyz_TurKinEneNonZero**0.5 &
!!$ & )
!!$ xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_StatStab &
!!$ & * xyz_ShHat &
!!$ & * xyz_TurKinEneNonZero**1.5
!!$ xyz_CBuo2 = - 2.0d0**1.5 * xyz_MixLength * xyz_StatStab &
!!$ & * ( xyz_DShHatDTKE &
!!$ & * xyz_TurKinEneNonZero**1.5 &
!!$ & + 1.5_DP * xyz_ShHat &
!!$ & * xyz_TurKinEneNonZero**0.5 &
!!$ & )
!!$ xyz_CDis1 = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$ & * xyz_TurKinEneNonZero**1.5
!!$ xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$ & * xyz_TurKinEneNonZero**0.5
xyz_CShe1 = sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq * xyz_Sm * sqrt( xyz_TurKinEneNonZero )
!!$ xyz_CShe2 = 0.5_DP * sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq &
!!$ & * xyz_Sm &
!!$ & / sqrt( xyz_TurKinEneNonZero )
xyz_CShe2 = 0.0_DP
xyz_CBuo1 = - sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab * xyz_Sh * sqrt( xyz_TurKinEneNonZero )
!!$ xyz_CBuo2 = - 0.5_DP * sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab &
!!$ & * xyz_Sh &
!!$ & / sqrt( xyz_TurKinEneNonZero )
xyz_CBuo2 = 0.0_DP
xyz_CDis1 = 2.0_DP**1.5_DP / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**1.5_DP
xyz_CDis2 = 1.5_DP * 2.0_DP**1.5_DP / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**0.5_DP
nloop = kmax
loop_solve : do iloop = 1, nloop
!
! Construct implicit matrix from transfer coefficient of vertical
! diffusion scheme (turbulent kinetic energy)
!
k = 1
xyza_TurKinEneMtx(:,:,k,-1) = 0.0_DP
xyza_TurKinEneMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_TurKinEneTransCoef(:,:,k-1) + xyr_TurKinEneTransCoef(:,:,k ) + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) )
xyza_TurKinEneMtx(:,:,k, 1) = - xyr_TurKinEneTransCoef(:,:,k )
!
do k = 2, kmax-1
xyza_TurKinEneMtx(:,:,k,-1) = - xyr_TurKinEneTransCoef(:,:,k-1)
xyza_TurKinEneMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_TurKinEneTransCoef(:,:,k-1) + xyr_TurKinEneTransCoef(:,:,k ) + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) )
xyza_TurKinEneMtx(:,:,k, 1) = - xyr_TurKinEneTransCoef(:,:,k )
end do
!
k = kmax
xyza_TurKinEneMtx(:,:,k,-1) = - xyr_TurKinEneTransCoef(:,:,k-1)
xyza_TurKinEneMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_TurKinEneTransCoef(:,:,k-1) + xyr_TurKinEneTransCoef(:,:,k ) + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) )
xyza_TurKinEneMtx(:,:,k, 1) = 0.0_DP
do k = 1, kmax
xyz_TurKinEneVec(:,:,k) = - ( xyr_TurKinEneFlux(:,:,k) - xyr_TurKinEneFlux(:,:,k-1) ) - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav * ( xyz_CShe1(:,:,k) + xyz_CBuo1(:,:,k) - xyz_CDis1(:,:,k) )
end do
!
! Solve simultaneous linear equations by use of LU decomposition technique
!
xyza_TurKinEneLUMtx = xyza_TurKinEneMtx
!
call PhyImplLUDecomp3( xyza_TurKinEneLUMtx, imax * jmax, kmax )
xyz_DelTurKinEneLUVec = xyz_TurKinEneVec
!
call PhyImplLUSolve3( xyz_DelTurKinEneLUVec, xyza_TurKinEneLUMtx, 1, imax * jmax , kmax )
xyz_DTurKinEneDt = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime )
! Calculation of dissipation rate of turbulent kinetic energy
!
! Calculate production rate of turbulent kinetic energy
! by shear and buoyancy
xyz_TurKinEneProShear = xyz_CShe1 + xyz_CShe2 * xyz_DTurKinEneDt * 2.0_DP * DelTime
xyz_TurKinEneProBuoya = xyz_CBuo1 + xyz_CBuo2 * xyz_DTurKinEneDt * 2.0_DP * DelTime
xyz_TurKinEneDiss = xyz_CDis1 + xyz_CDis2 * xyz_DTurKinEneDt * 2.0_DP * DelTime
! Check of turbulent kinetic energy dissipation rate
! If it is negative, tendency is recalculated without dissipation.
!
FlagReCalc = .false.
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_TurKinEneDiss(i,j,k) < 0.0_DP ) then
xyz_CDis1(i,j,k) = 0.0_DP
xyz_CDis2(i,j,k) = 0.0_DP
FlagReCalc = .true.
end if
end do
end do
end do
! Check convergence
a_FlagReCalcLocal = FlagReCalc
call MPIWrapperChkTrue( 1, a_FlagReCalcLocal, a_FlagReCalcGlobal )
if ( .not. a_FlagReCalcGlobal(1) ) exit loop_solve
end do loop_solve
!!$ write( 6, * ) TimeN, iloop
!!$ write( 6, * ) xyz_TurKinEne(0,jmax/2+1,1:4)
!!$ write( 6, * ) xyz_TempDiffCoef(0,jmax/2+1,1:4)
!!$ write( 6, * ) xyz_TurKinEneProShear(0,jmax/2+1,1:4)
!!$ write( 6, * ) xyz_TurKinEneProBuoya(0,jmax/2+1,1:4)
!!$ write( 6, * ) xyz_TurKinEneDiss(0,jmax/2+1,1:4)
!!$ write( 6, * ) xyz_DTurKinEneDt(0,jmax/2+1,1:4)
! 拡散係数の出力
! Output diffusion coefficients
!
! 拡散係数出力
! Diffusion coeffficients output
!
call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef )
call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef )
call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef )
call HistoryAutoPut( TimeN, 'TKEPShear', xyz_TurKinEneProShear )
call HistoryAutoPut( TimeN, 'TKEPBuoy' , xyz_TurKinEneProBuoya )
call HistoryAutoPut( TimeN, 'TKEDiss' , xyz_TurKinEneDiss )
call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine VDiffusionMY25
| Subroutine : | |||
| 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)
| ||
| xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | 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)
| ||
| 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)
| ||
| xyr_Height(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)
| ||
| xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_SurfMomFluxX(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfMomFluxY(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyz_DTurKinEneDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated by use of MY2.5 model.
subroutine VDiffusionMY251DWrapper3D( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyz_TurKinEne, xy_SurfMomFluxX, xy_SurfMomFluxY, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyz_DTurKinEneDt )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated by use of MY2.5 model.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm, Grav, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! $ u $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! $ v $ . 南北風速. Northward wind
real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q $ . 質量混合比. Mass mixing ratio
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T} $ . 温度 (半整数レベル).
! Temperature (half level)
real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
! $ T_v $ . 仮温度. Virtual temperature
real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air 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):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
! 高度 (半整数レベル).
! Height (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):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax)
!
! Turbulent kinetic energy (m2 s-2)
real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax)
!
! Eastward momentum flux at surface
real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax)
!
! Northward momentum flux at surface
real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 質量フラックス.
! Mass flux of compositions
real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Diffusion coefficient: temperature
real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax)
!
! Tendency of turbulent kinetic energy
! 作業変数
! Work variables
!
real(DP) :: z_U (1:kmax)
real(DP) :: z_V (1:kmax)
real(DP) :: zf_QMix(1:kmax, 1:ncmax)
real(DP) :: z_Temp (1:kmax)
real(DP) :: r_Temp (0:kmax)
real(DP) :: z_VirTemp (1:kmax)
real(DP) :: r_VirTemp (0:kmax)
real(DP) :: r_Press (0:kmax)
real(DP) :: SurfHeight
real(DP) :: z_Height (1:kmax)
real(DP) :: r_Height (0:kmax)
real(DP) :: z_Exner (1:kmax)
real(DP) :: r_Exner (0:kmax)
real(DP) :: z_TurKinEne(1:kmax)
real(DP) :: SurfMomFluxX
real(DP) :: SurfMomFluxY
real(DP) :: r_MomFluxX (0:kmax)
real(DP) :: r_MomFluxY (0:kmax)
real(DP) :: r_HeatFlux (0:kmax)
real(DP) :: rf_QMixFlux(0:kmax, 1:ncmax)
real(DP) :: r_VelDiffCoef (0:kmax)
real(DP) :: r_TempDiffCoef(0:kmax)
real(DP) :: r_QMixDiffCoef(0:kmax)
real(DP) :: z_DTurKinEneDt (1:kmax)
!
! Tendency of turbulent kinetic energy
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
do j = 1, jmax
do i = 0, imax-1
z_U = xyz_U (i,j,:)
z_V = xyz_V (i,j,:)
zf_QMix = xyzf_QMix (i,j,:,:)
z_Temp = xyz_Temp (i,j,:)
r_Temp = xyr_Temp (i,j,:)
z_VirTemp = xyz_VirTemp (i,j,:)
r_VirTemp = xyr_VirTemp (i,j,:)
r_Press = xyr_Press (i,j,:)
SurfHeight = xy_SurfHeight (i,j)
z_Height = xyz_Height (i,j,:)
r_Height = xyr_Height (i,j,:)
z_Exner = xyz_Exner (i,j,:)
r_Exner = xyr_Exner (i,j,:)
z_TurKinEne = xyz_TurKinEne (i,j,:)
SurfMomFluxX = xy_SurfMomFluxX(i,j)
SurfMomFluxY = xy_SurfMomFluxY(i,j)
call VDiffusionMY251D( z_U, z_V, zf_QMix, z_Temp, r_Temp, z_VirTemp, r_VirTemp, r_Press, SurfHeight, z_Height, r_Height, z_Exner, r_Exner, z_TurKinEne, SurfMomFluxX, SurfMomFluxY, r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux, r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, z_DTurKinEneDt )
xyr_MomFluxX (i,j,:) = r_MomFluxX
xyr_MomFluxY (i,j,:) = r_MomFluxY
xyr_HeatFlux (i,j,:) = r_HeatFlux
xyrf_QMixFlux (i,j,:,:) = rf_QMixFlux
xyr_VelDiffCoef (i,j,:) = r_VelDiffCoef
xyr_TempDiffCoef(i,j,:) = r_TempDiffCoef
xyr_QMixDiffCoef(i,j,:) = r_QMixDiffCoef
xyz_DTurKinEneDt (i,j,:) = z_DTurKinEneDt
end do
end do
!!$ ! 拡散係数の出力
!!$ ! Output diffusion coefficients
!!$ !
!!$
!!$ ! 拡散係数出力
!!$ ! Diffusion coeffficients output
!!$ !
!!$ call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef )
!!$ call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef )
!!$ call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef )
!!$
!!$ call HistoryAutoPut( TimeN, 'TKEPShear', xyz_TurKinEneProShear )
!!$ call HistoryAutoPut( TimeN, 'TKEPBuoy' , xyz_TurKinEneProBuoya )
!!$ call HistoryAutoPut( TimeN, 'TKEDiss' , xyz_TurKinEneDiss )
!!$
!!$ call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine VDiffusionMY251DWrapper3D
| Subroutine : | |||
| 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)
| ||
| xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | 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)
| ||
| 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)
| ||
| xyr_Height(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)
| ||
| xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_SurfMomFluxX(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfMomFluxY(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyz_DTurKinEneDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated by use of MY2.5 model.
subroutine VDiffusionMY25GBT94( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyz_TurKinEne, xy_SurfMomFluxX, xy_SurfMomFluxY, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyz_DTurKinEneDt )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated by use of MY2.5 model.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm, Grav, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! MPI 関連ルーチン
! MPI related routines
!
use mpi_wrapper, only : MPIWrapperChkTrue
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! $ u $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! $ v $ . 南北風速. Northward wind
real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q $ . 質量混合比. Mass mixing ratio
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T} $ . 温度 (半整数レベル).
! Temperature (half level)
real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
! $ T_v $ . 仮温度. Virtual temperature
real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air 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):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
! 高度 (半整数レベル).
! Height (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):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax)
!
! Turbulent kinetic energy (m2 s-2)
real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax)
!
! Eastward momentum flux at surface
real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax)
!
! Northward momentum flux at surface
real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 質量フラックス.
! Mass flux of compositions
real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Diffusion coefficient: temperature
real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax)
!
! Tendency of turbulent kinetic energy
! 作業変数
! Work variables
!
real(DP) :: xyz_MixLength(0:imax-1, 1:jmax, 1:kmax)
! 混合距離.
! Mixing length
real(DP) :: xyz_DVelDzSq(0:imax-1, 1:jmax, 1:kmax)
!
! Vertical shear squared (s-2)
real(DP) :: xyz_StatStab(0:imax-1, 1:jmax, 1:kmax)
!
! Static stability (s-2)
!!$ real(DP) :: GhMin
!!$ !
!!$ ! Minimum of G_h
!!$ real(DP) :: GhMax
!!$ !
!!$ ! Maximum of G_h
real(DP) :: xyz_Gm(0:imax-1, 1:jmax, 1:kmax)
!
! G_m
real(DP) :: xyz_Gh(0:imax-1, 1:jmax, 1:kmax)
!
! G_h
!!$ real(DP) :: xyz_SmHat(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! \hat{S}_M
!!$ real(DP) :: xyz_ShHat(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! \hat{S}_h
!!$ real(DP) :: xyz_DSmHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! derivative of \hat{S}_M
!!$ real(DP) :: xyz_DShHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$ !
!!$ ! derivative of \hat{S}_h
real(DP) :: xyz_Sm(0:imax-1, 1:jmax, 1:kmax)
!
! S_M
real(DP) :: xyz_Sh(0:imax-1, 1:jmax, 1:kmax)
!
! S_h
real(DP), parameter :: Stke = 0.2_DP
!
! S_{TKE} = 0.2
real(DP) :: xyz_VelDiffCoef (0:imax-1, 1:jmax, 1:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP) :: xyz_TempDiffCoef (0:imax-1, 1:jmax, 1:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_TurKinEneDiffCoef (0:imax-1, 1:jmax, 0:kmax)
!
! Diffusion coefficient: turbulent kinetic energy
real(DP) :: xyz_TurKinEneDiffCoef (0:imax-1, 1:jmax, 1:kmax)
!
! Diffusion coefficient: turbulent kinetic energy
real(DP) :: xyr_TurKinEneTransCoef(0:imax-1, 1:jmax, 0:kmax)
!
! Transfer coefficient: turbulent kinetic energy
real(DP) :: xyr_TurKinEneFlux(0:imax-1, 1:jmax, 0:kmax)
!
! Turbulent energy flux
real(DP) :: xy_FricVelSq (0:imax-1, 1:jmax)
real(DP) :: xy_TurKinEneAtLB(0:imax-1, 1:jmax)
real(DP) :: xyza_TurKinEneMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
!
! Implicit matrix for turbulent kinetic energy
real(DP) :: xyz_TurKinEneVec(0:imax-1, 1:jmax, 1:kmax)
!
! Implicit vector for turbulent kinetic energy
real(DP) :: xyza_TurKinEneLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP) :: xyz_DelTurKinEneLUVec(0:imax-1, 1:jmax, 1:kmax)
!
! Tendency of turbulent kinetic energy
real(DP) :: xyz_TurKinEneDiss(0:imax-1, 1:jmax, 1:kmax)
!
! Dissipation rate of turbulent kinetic energy (m2 s-3)
real(DP) :: xyz_TurKinEneNonZero(0:imax-1, 1:jmax, 1:kmax)
!
! Turbulent kinetic energy with offset (m2 s-2)
real(DP), parameter :: TurKinEneOffset = ( 1.0e-3_DP )**2 / 2.0_DP
logical :: xyz_FlagUseRiNum (0:imax-1, 1:jmax, 1:kmax)
logical :: xyz_FlagTKEAsymptToZero(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_RiNum(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: Beta1
real(DP) :: Beta2
real(DP) :: Beta3
real(DP) :: Beta4
real(DP) :: MixLength
real(DP) :: xyz_A1 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_A2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_R1 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_R2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_SqrtA1(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_TKEInit (0:imax-1, 1:jmax, 1:kmax)
!
! Turbulent kinetic energy for current time step
! with offset (m2 s-2)
real(DP) :: xyz_TKEx2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_SqrtTKEx2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_TKETentative(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: Alpha
real(DP) :: BetaSq
real(DP) :: StatStab
real(DP) :: RiNum
real(DP) :: DVelDzSq
real(DP), parameter :: Epsilon = 1.0e-10_DP
real(DP), parameter :: CrtlRiNum = 0.195_DP
real(DP), parameter :: CrtlShear = 0.001_DP / 1000.0_DP
!!$ logical :: FlagReCalc
!!$ !
!!$ ! Flag for recalculation
!!$ logical :: a_FlagReCalcLocal (1)
!!$ logical :: a_FlagReCalcGlobal(1)
!!$ integer :: iloop
!!$ integer :: nloop
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
!!$ integer:: n ! 組成方向に回る DO ループ用作業変数
!!$ ! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
!
! Calculation of vertical shear squared
do k = 1, kmax
if ( k == 1 ) then
xyz_DVelDzSq(:,:,k) = ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k ) )**2 + ( xyz_V(:,:,k+1) - xyz_V(:,:,k ) )**2 ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) )**2
!!$ xyz_DVelDzSq(:,:,k) = &
!!$ & ( ( xyz_U(:,:,k+1) - 0.0_DP )**2 &
!!$ & + ( xyz_V(:,:,k+1) - 0.0_DP )**2 ) &
!!$ & / ( xyz_Height(:,:,k+1) - xy_SurfHeight )**2
else if ( k == kmax ) then
xyz_DVelDzSq(:,:,k) = ( ( xyz_U(:,:,k ) - xyz_U(:,:,k-1) )**2 + ( xyz_V(:,:,k ) - xyz_V(:,:,k-1) )**2 ) / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) )**2
else
xyz_DVelDzSq(:,:,k) = ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k-1) )**2 + ( xyz_V(:,:,k+1) - xyz_V(:,:,k-1) )**2 ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )**2
end if
end do
! Calculation of static stability
do k = 1, kmax
if ( k == 1 ) then
xyz_StatStab(:,:,k) = Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) )
else if ( k == kmax ) then
xyz_StatStab(:,:,k) = Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) * ( xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) )
else
xyz_StatStab(:,:,k) = Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )
end if
end do
! 混合距離の算出
! Calculate mixing length
!
do k = 1, kmax
xyz_MixLength(:,:,k) = FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / (1.0_DP + FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / MixLengthMax )
end do
!!$ ! Limit mixing length (Galperin et al., 1988) and avoid zero
!!$ xyz_MixLength = &
!!$ & min( xyz_MixLength, &
!!$ & 0.53_DP &
!!$ & * sqrt( 2.0_DP * xyz_TurKinEneNonZero / max( xyz_StatStab, 1.0d-10 ) ) ) &
!!$ & + 1.0d-10
!********************************************************************
! Gerrity et al. (1994)
! Set flag for using Richardson number
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_DVelDzSq(i,j,k) < CrtlShear**2 ) then
xyz_FlagUseRiNum(i,j,k) = .false.
else
xyz_FlagUseRiNum(i,j,k) = .true.
end if
end do
end do
end do
! Calculation of Richardson number
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FlagUseRiNum(i,j,k) ) then
xyz_RiNum(i,j,k) = xyz_StatStab(i,j,k) / xyz_DVelDzSq(i,j,k)
else
xyz_RiNum(i,j,k) = - 1.0e100_DP
end if
end do
end do
end do
! Set flag for selecting an asymptotic equation
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FlagUseRiNum(i,j,k) ) then
if ( CrtlRiNum <= xyz_RiNum(i,j,k) ) then
! Ric <= Ri
xyz_FlagTKEAsymptToZero = .true.
else
! Ri < Ric
xyz_FlagTKEAsymptToZero = .false.
end if
else
if ( xyz_StatStab(i,j,k) > 0.0_DP ) then
xyz_FlagTKEAsymptToZero = .true.
else
xyz_FlagTKEAsymptToZero = .false.
end if
end if
end do
end do
end do
! Calculation of roos for equations for steady state
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FlagUseRiNum(i,j,k) ) then
! Calculation with Richardson number
MixLength = xyz_MixLength(i,j,k)
DVelDzSq = xyz_DVelDzSq (i,j,k)
RiNum = xyz_RiNum (i,j,k)
!
Beta1 = - MixLength**2 * DVelDzSq * ( 6.53_DP - 49.0_DP * RiNum )
Beta2 = - MixLength**4 * DVelDzSq**2 * ( 51.2_DP - 262.7_DP * RiNum ) * RiNum
Beta3 = MixLength**2 * DVelDzSq * ( 5.08_DP + 36.7_DP * RiNum )
Beta4 = MixLength**4 * DVelDzSq**2 * ( 88.8_DP + 187.4_DP * RiNum ) * RiNum
!
xyz_A1(i,j,k) = - Beta1 / 2.0_DP + sqrt( Beta1**2 - 4.0_DP * Beta2 ) / 2.0_DP
xyz_A2(i,j,k) = - Beta1 / 2.0_DP - sqrt( Beta1**2 - 4.0_DP * Beta2 ) / 2.0_DP
xyz_R1(i,j,k) = - Beta3 / 2.0_DP + sqrt( Beta3**2 - 4.0_DP * Beta4 ) / 2.0_DP
xyz_R2(i,j,k) = - Beta3 / 2.0_DP - sqrt( Beta3**2 - 4.0_DP * Beta4 ) / 2.0_DP
else
! Calculation without Richardson number
MixLength = xyz_MixLength(i,j,k)
StatStab = xyz_StatStab (i,j,k)
!
xyz_A1(i,j,k) = 24.49_DP * MixLength**2 * ( - StatStab ) + 18.36_DP * MixLength**2 * abs( - StatStab )
xyz_A2(i,j,k) = 24.49_DP * MixLength**2 * ( - StatStab ) - 18.36_DP * MixLength**2 * abs( - StatStab )
xyz_R1(i,j,k) = 18.35_DP * MixLength**2 * ( - StatStab ) + 12.22_DP * MixLength**2 * abs( - StatStab )
xyz_R2(i,j,k) = 18.35_DP * MixLength**2 * ( - StatStab ) - 12.22_DP * MixLength**2 * abs( - StatStab )
end if
end do
end do
end do
! Set turbulent kinetic energy at current time step
!
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FlagUseRiNum(i,j,k) ) then
if ( xyz_RiNum(i,j,k) < 0.0_DP ) then
! Ri < 0 (or CT > 0)
xyz_TKEInit(i,j,k) = max( xyz_TurKinEne(i,j,k), xyz_R1(i,j,k) / 2.0_DP )
else if ( xyz_RiNum(i,j,k) < CrtlRiNum ) then
! 0 <= Ri < Ric
xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k)
else
! Ric <= Ri
xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k)
end if
else
if ( xyz_StatStab(i,j,k) < 0.0_DP ) then
! CT > 0
xyz_TKEInit(i,j,k) = max( xyz_TurKinEne(i,j,k), xyz_R1(i,j,k) / 2.0_DP )
else
! CT <= 0
xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k)
end if
end if
xyz_TKEInit(i,j,k) = xyz_TKEInit(i,j,k) + TurKinEneOffset
end do
end do
end do
!
xyz_TKEx2 = 2.0_DP * xyz_TKEInit
xyz_SqrtTKEx2 = sqrt( xyz_TKEx2 )
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FlagTKEAsymptToZero(i,j,k) ) then
xyz_SqrtA1(i,j,k) = 1.0e100_DP
else
xyz_SqrtA1(i,j,k) = sqrt( xyz_A1(i,j,k) )
end if
end do
end do
end do
!
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FlagTKEAsymptToZero(i,j,k) ) then
! Use equations asymptotic to zero TKE, Eq. (10)
! Eq. (11)
BetaSq = ( xyz_TKEx2(i,j,k) - xyz_A1(i,j,k)**2 ) * ( xyz_TKEx2(i,j,k) - xyz_A2(i,j,k)**2 ) / ( MYConstB1 * xyz_MixLength(i,j,k) * ( xyz_TKEx2(i,j,k) - xyz_R1(i,j,k) ) * ( xyz_TKEx2(i,j,k) - xyz_R2(i,j,k) ) )
! Eq. (10)
xyz_TKETentative(i,j,k) = xyz_TKEInit(i,j,k) / 2.0_DP / ( 1.0_DP + xyz_SqrtTKEx2(i,j,k) * BetaSq * ( 2.0_DP * DelTime ) )**2
else
! Use equations asymptotic to certain TKE, Eq. (8)
! Eq. (7b)
BetaSq = xyz_TKEx2(i,j,k) * ( xyz_TKEx2(i,j,k) - xyz_A2(i,j,k)**2 ) / ( MYConstB1 * xyz_MixLength(i,j,k) * ( xyz_TKEx2(i,j,k) - xyz_R1(i,j,k) ) * ( xyz_TKEx2(i,j,k) - xyz_R2(i,j,k) ) )
! Eq. (9)
Alpha = ( xyz_SqrtTKEx2(i,j,k) - xyz_SqrtA1(i,j,k) ) / ( xyz_SqrtTKEx2(i,j,k) + xyz_SqrtA1(i,j,k) ) * exp( - 2.0_DP * xyz_SqrtA1(i,j,k) * BetaSq * ( 2.0_DP * DelTime ) )
Alpha = max( min( Alpha, 1.0_DP - Epsilon ), -1.0_DP )
! Eq. (8)
xyz_TKETentative(i,j,k) = xyz_A1(i,j,k) / 2.0_DP * ( ( 1.0_DP + Alpha ) / ( 1.0_DP - Alpha ) )**2
end if
end do
end do
end do
!********************************************************************
! Set turbulent kinetic energy for diffusion calculation
!
!!$ xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset
xyz_TurKinEneNonZero = xyz_TKETentative + TurKinEneOffset
xyz_Gh = - xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_StatStab
! Actually, xyz_Gm is not used below.
xyz_Gm = xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_DVelDzSq
! Limit Gh (Galperin et al., 1988)
!!$ GhMin = - 0.53d0**2
!!$ GhMax = 1.0_DP &
!!$ & / ( MYConstA2 &
!!$ & * ( 12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) )
!!$ xyz_Gh = max( GhMin, min( xyz_Gh, GhMax ) )
xyz_Sh = MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) / ( 1.0_DP - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_Gh )
xyz_Sm = ( MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 - 6.0_DP * MYConstA1 / MYConstB1 ) + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) * xyz_Gh * xyz_Sh ) / ( 1.0_DP - 9.0_DP * MYConstA1 * MYConstA2 * xyz_Gh )
! 拡散係数の計算
! Calculation of diffusion coefficient
!
xyz_VelDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm
xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh
!
do k = 0, kmax
if ( ( k == 0 ) .or. ( k == kmax ) ) then
xyr_VelDiffCoef (:,:,k) = 0.0_DP
xyr_TempDiffCoef(:,:,k) = 0.0_DP
else
xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP
xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP
end if
end do
!
do k = 1, kmax-1
do j = 1, jmax
do i = 0, imax-1
xyr_VelDiffCoef(i,j,k) = max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), VelDiffCoefMin )
xyr_TempDiffCoef(i,j,k) = max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), TempDiffCoefMin )
end do
end do
end do
!
xyr_QMixDiffCoef = xyr_TempDiffCoef
! 輸送係数とフラックスの計算
! Calculate transfer coefficient and flux
!
call VDiffusionCalcFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux )
! Calculate tendency of turbulent kinetic energy
! Set diffusion coefficient for turbulent kinetic energy
xyz_TurKinEneDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * Stke
!
do k = 0, kmax
if ( k == 0 ) then
xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1)
else if ( k == kmax ) then
xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,kmax)
else
xyr_TurKinEneDiffCoef(:,:,k) = ( xyz_TurKinEneDiffCoef(:,:,k) + xyz_TurKinEneDiffCoef(:,:,k+1) ) / 2.0_DP
end if
end do
! Calculate turbulent kinetic energy at lower boundary
!
xy_FricVelSq = sqrt( xy_SurfMomFluxX**2 + xy_SurfMomFluxY**2 ) / ( xyr_Press(:,:,0) / ( GasRDry * xyr_VirTemp(:,:,0) ) )
xy_TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * xy_FricVelSq
xy_TurKinEneAtLB = xy_TurKinEneAtLB + TurKinEneOffset
! Calculate transfer coefficient and flux of turbulent kinetic energy
!
! When transfer coefficient at lower boundary is calculated,
! diffusion coefficient at mid-point of 1st layer is used.
! In addition, transfer coefficient at upper boundary is assumed
! to be zero.
k = 0
xyr_TurKinEneTransCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,1) - xy_SurfHeight )
do k = 1, kmax-1
xyr_TurKinEneTransCoef(:,:,k) = xyr_TurKinEneDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
end do
k = kmax
xyr_TurKinEneTransCoef(:,:,k) = 0.0_DP
!
do k = 1, kmax-1
xyr_TurKinEneFlux(:,:,k) = - xyr_TurKinEneTransCoef(:,:,k) * ( xyz_TurKinEneNonZero(:,:,k+1) - xyz_TurKinEneNonZero(:,:,k) )
end do
k = 0
xyr_TurKinEneFlux(:,:,k) = - xyr_TurKinEneTransCoef(:,:,k) * ( xyz_TurKinEneNonZero(:,:,k+1) - xy_TurKinEneAtLB )
k = kmax
xyr_TurKinEneFlux(:,:,k) = 0.0_DP
!
! Construct implicit matrix from transfer coefficient of vertical
! diffusion scheme (turbulent kinetic energy)
!
k = 1
xyza_TurKinEneMtx(:,:,k,-1) = 0.0_DP
xyza_TurKinEneMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_TurKinEneTransCoef(:,:,k-1) + xyr_TurKinEneTransCoef(:,:,k )
xyza_TurKinEneMtx(:,:,k, 1) = - xyr_TurKinEneTransCoef(:,:,k )
!
do k = 2, kmax-1
xyza_TurKinEneMtx(:,:,k,-1) = - xyr_TurKinEneTransCoef(:,:,k-1)
xyza_TurKinEneMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_TurKinEneTransCoef(:,:,k-1) + xyr_TurKinEneTransCoef(:,:,k )
xyza_TurKinEneMtx(:,:,k, 1) = - xyr_TurKinEneTransCoef(:,:,k )
end do
!
k = kmax
xyza_TurKinEneMtx(:,:,k,-1) = - xyr_TurKinEneTransCoef(:,:,k-1)
xyza_TurKinEneMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_TurKinEneTransCoef(:,:,k-1) + xyr_TurKinEneTransCoef(:,:,k )
xyza_TurKinEneMtx(:,:,k, 1) = 0.0_DP
do k = 1, kmax
xyz_TurKinEneVec(:,:,k) = - ( xyr_TurKinEneFlux(:,:,k) - xyr_TurKinEneFlux(:,:,k-1) )
end do
!
! Solve simultaneous linear equations by use of LU decomposition technique
!
xyza_TurKinEneLUMtx = xyza_TurKinEneMtx
!
call PhyImplLUDecomp3( xyza_TurKinEneLUMtx, imax * jmax, kmax )
xyz_DelTurKinEneLUVec = xyz_TurKinEneVec
!
call PhyImplLUSolve3( xyz_DelTurKinEneLUVec, xyza_TurKinEneLUMtx, 1, imax * jmax , kmax )
!!$ xyz_DTurKinEneDt = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime )
xyz_TKETentative = xyz_TKETentative + xyz_DelTurKinEneLUVec
xyz_DTurKinEneDt = ( xyz_TKETentative - xyz_TurKinEne ) / ( 2.0_DP * DelTime )
!!$ write( 6, * ) TimeN, iloop
!!$ write( 6, * ) xyz_TurKinEne(0,jmax/2+1,1:4)
!!$ write( 6, * ) xyz_TempDiffCoef(0,jmax/2+1,1:4)
!!$ write( 6, * ) xyz_DTurKinEneDt(0,jmax/2+1,1:4)
! 拡散係数の出力
! Output diffusion coefficients
!
! 拡散係数出力
! Diffusion coeffficients output
!
call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef )
call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef )
call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef )
call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine VDiffusionMY25GBT94
| Subroutine : | |||
| 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)
| ||
| 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)
| ||
| 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)
| ||
| xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_Height(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | 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 VDiffusionOutPut( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )
!
! フラックス (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 constants settings
!
use constants, only: CpDry, LatentHeat, GasRDry
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
! 時刻管理
! 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 wind 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)
! 質量フラックス.
! Mass flux of constituents
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} $ . 混合比変化.
! Mass mixing ratio 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):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
! 高度 (整数レベル).
! Height (full level)
real(DP), intent(in):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
! 出力のための作業変数
! 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)
! 質量フラックス.
! Mass flux of constituents
! 作業変数
! Work variables
!
real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP) :: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass of constituents
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
real(DP):: LCp
! $ L / C_p $ [K].
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_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 transfer coefficient
!
xyr_VelTransCoef (:,:,0) = 0.0_DP
xyr_VelTransCoef (:,:,kmax) = 0.0_DP
xyr_TempTransCoef(:,:,0) = 0.0_DP
xyr_TempTransCoef(:,:,kmax) = 0.0_DP
xyr_QMixTransCoef(:,:,0) = 0.0_DP
xyr_QMixTransCoef(:,:,kmax) = 0.0_DP
do k = 1, kmax-1
xyr_VelTransCoef(:,:,k) = xyr_VelDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
xyr_TempTransCoef(:,:,k) = xyr_TempDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
xyr_QMixTransCoef(:,:,k) = xyr_QMixDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
end do
! 風速, 温度, 比湿フラックス補正
! Correct fluxes of wind, temperature, specific humidity
!
LCp = LatentHeat / CpDry
do k = 1, kmax-1
do j = 1, jmax
do i = 0, imax-1
xyr_MomFluxXCor( i,j,k ) = xyr_MomFluxX( i,j,k ) + ( xyz_DUDt( i,j,k ) - xyz_DUDt( i,j,k+1 ) ) * xyr_VelTransCoef( i,j,k ) * DelTime
xyr_MomFluxYCor( i,j,k ) = xyr_MomFluxY( i,j,k ) + ( xyz_DVDt( i,j,k ) - xyz_DVDt( i,j,k+1 ) ) * xyr_VelTransCoef( i,j,k ) * DelTime
xyr_HeatFluxCor( i,j,k ) = xyr_HeatFlux( i,j,k ) + ( xyz_DTempDt( i,j,k ) / xyz_Exner( i,j,k ) - xyz_DTempDt( i,j,k+1 ) / xyz_Exner( i,j,k+1 ) ) * CpDry * xyr_TempTransCoef( i,j,k ) * xyr_Exner( i,j,k ) * DelTime
end do
end do
end do
do n = 1, ncmax
do k = 1, kmax-1
do j = 1, jmax
do i = 0, imax-1
xyrf_QMixFluxCor( i,j,k,n ) = xyrf_QMixFlux( i,j,k,n ) + ( xyzf_DQMixDt( i,j,k,n ) - xyzf_DQMixDt( i,j,k+1,n ) ) * CpDry * xyr_QMixTransCoef( i,j,k ) * LCp * DelTime
end do
end do
end do
end do
xyr_MomFluxXCor (:,:,0) = 0.0_DP
xyr_MomFluxXCor (:,:,kmax) = 0.0_DP
xyr_MomFluxYCor (:,:,0) = 0.0_DP
xyr_MomFluxYCor (:,:,kmax) = 0.0_DP
xyr_HeatFluxCor (:,:,0) = 0.0_DP
xyr_HeatFluxCor (:,:,kmax) = 0.0_DP
do n = 1, ncmax
xyrf_QMixFluxCor(:,:,0,n) = 0.0_DP
xyrf_QMixFluxCor(:,:,kmax,n) = 0.0_DP
end do
! MEMO
! Output values of surface fluxes in MomFluxX, MomFluxY, HeatFlux, and QVapFlux
! are not correct. (YOT, 2009/08/14)
! Please refer to output variables, 'TauX', 'TauY', 'Sens', and 'Evap' for those
! values. (YOT, 2011/05/28)
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'MomFluxX', xyr_MomFluxXCor )
call HistoryAutoPut( TimeN, 'MomFluxY', xyr_MomFluxYCor )
call HistoryAutoPut( TimeN, 'HeatFlux', xyr_HeatFluxCor )
call HistoryAutoPut( TimeN, 'QVapFlux', xyrf_QMixFluxCor )
call HistoryAutoPut( TimeN, 'DUDtVDiff' , xyz_DUDt )
call HistoryAutoPut( TimeN, 'DVDtVDiff' , xyz_DVDt )
call HistoryAutoPut( TimeN, 'DTempDtVDiff', xyz_DTempDt )
call HistoryAutoPut( TimeN, 'DQVapDtVDiff', xyzf_DQMixDt(:,:,:,IndexH2OVap) )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine VDiffusionOutPut
| Variable : | |||
| BulkRiNumMin : | real(DP), save
|
| Variable : | |||
| ConstDiffCoefH : | real(DP), save
|
| Variable : | |||
| ConstDiffCoefM : | real(DP), save
|
| Variable : | |||
| FlagConstDiffCoef : | logical , save
|
| Variable : | |||
| TempDiffCoefMax : | real(DP), save
|
| Variable : | |||
| TempDiffCoefMin : | real(DP), save
|
| Subroutine : | |||
| xy_SurfHeight(0:imax-1,1:jmax) : | real(DP), intent(in)
| ||
| xyr_Height(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_DVelDz(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_BulkRiNum(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated.
subroutine VDiffCoefficient( xy_SurfHeight, xyr_Height, xyr_DVelDz, xyr_BulkRiNum, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated.
!
! モジュール引用 ; USE statements
!
! 時刻管理
! Time control
!
use timeset, only: TimeN, TimesetClockStart, TimesetClockStop
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm
! $ k $ .
! カルマン定数.
! Karman constant
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xy_SurfHeight (0:imax-1,1:jmax)
! $ z_s $ . 地表面高度.
! Surface height.
real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
! 高度 (半整数レベル).
! Height (half level)
real(DP), intent(in):: xyr_DVelDz (0:imax-1, 1:jmax, 0:kmax)
! $ \DD{|\Dvect{v}|}{z} $
real(DP), intent(in):: xyr_BulkRiNum (0:imax-1, 1:jmax, 0:kmax)
! バルク $ R_i $ 数.
! Bulk $ R_i $
real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：質量
! Diffusion coefficient: mass of constituents
! 作業変数
! Work variables
!
real(DP):: xyr_FluxRiNum (0:imax-1, 1:jmax, 0:kmax)
! フラックス $ R_i $ 数.
! Flux $ R_i $ number
real(DP):: xyr_Sh (0:imax-1, 1:jmax, 0:kmax)
! $ S_h $ (温度, 比湿).
! $ S_h $ (temperature, specific humidity)
real(DP):: xyr_Sm (0:imax-1, 1:jmax, 0:kmax)
! $ S_m $ (運動量).
! $ S_m $ (momentum)
real(DP):: xyr_MixLength (0:imax-1, 1:jmax, 0:kmax)
! 混合距離.
! Mixing length
real(DP):: Alpha1, Alpha2
real(DP):: Beta1, Beta2, Beta3, Beta4
real(DP):: Gamma1, Gamma2
real(DP):: CrtlFluxRiNum
real(DP):: xyr_TurKinEne(0:imax-1, 1:jmax, 0:kmax)
!
! Turbulent kinetic energy
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 定数計算
! Calculate constants
!
Gamma1 = ( 1.0_DP / 3.0_DP ) - ( 2.0_DP * MYConstA1 / MYConstB1 )
Gamma2 = ( MYConstB2 / MYConstB1 ) + ( 6.0_DP * MYConstA1 / MYConstB1 )
Alpha1 = 3.0_DP * MYConstA2 * Gamma1
Alpha2 = 3.0_DP * MYConstA2 * ( Gamma1 + Gamma2 )
Beta1 = MYConstA1 * MYConstB1 * ( Gamma1 - MYConstC1 )
Beta2 = MYConstA1 * ( MYConstB1 * ( Gamma1 - MYConstC1 ) + 6.0_DP * MYConstA1 + 3.0_DP * MYConstA2 )
Beta3 = MYConstA2 * MYConstB1 * Gamma1
Beta4 = MYConstA2 * ( MYConstB1 * ( Gamma1 + Gamma2 ) - 3.0_DP * MYConstA1 )
CrtlFluxRiNum = Gamma1 / ( Gamma1 + Gamma2 )
! フラックス $ R_i $ 数の算出
! Calculate flux $ R_i $ number
!
xyr_FluxRiNum = ( Beta1 + Beta4 * xyr_BulkRiNum - sqrt( ( Beta1 + Beta4 * xyr_BulkRiNum )**2 - 4.0_DP * Beta2 * Beta3 * xyr_BulkRiNum ) ) / ( 2.0_DP * Beta2 )
! $ \tilde{S_h} $ と $ \tilde{S_m} $ の算出
! Calculate $ \tilde{S_h} $ and $ \tilde{S_m} $
!
xyr_Sh(:,:,kmax) = 0.0_DP
xyr_Sm(:,:,kmax) = 0.0_DP
do k = 0, kmax-1
do i = 0, imax-1
do j = 1, jmax
if ( xyr_FluxRiNum(i,j,k) < CrtlFluxRiNum ) then
xyr_Sh(i,j,k) = ( Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) / ( 1.0_DP - 1.0_DP * xyr_FluxRiNum(i,j,k) )
xyr_Sm(i,j,k) = ( Beta1 - Beta2 * xyr_FluxRiNum(i,j,k) ) / ( Beta3 - Beta4 * xyr_FluxRiNum(i,j,k) ) * xyr_Sh(i,j,k)
xyr_Sh(i,j,k) = max( xyr_Sh(i,j,k), ShMin )
xyr_Sm(i,j,k) = max( xyr_Sm(i,j,k), SmMin )
else
xyr_Sh(i,j,k) = ShMin
xyr_Sm(i,j,k) = SmMin
end if
end do
end do
end do
! 混合距離の算出
! Calculate mixing length
!
do k = 0, kmax
xyr_MixLength(:,:,k) = FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / (1.0_DP + FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / MixLengthMax )
end do
! 拡散係数の算出
! Calculate diffusion constants
!
xyr_VelDiffCoef = xyr_MixLength**2 * xyr_DVelDz * sqrt ( MYConstB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm ) * xyr_Sm
xyr_TempDiffCoef = xyr_MixLength ** 2 * xyr_DVelDz * sqrt ( MYConstB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm ) * xyr_Sh
do k = 0, kmax-1
do i = 0, imax-1
do j = 1, jmax
xyr_VelDiffCoef(i,j,k) = max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), VelDiffCoefMin )
xyr_TempDiffCoef(i,j,k) = max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), TempDiffCoefMin )
end do
end do
end do
xyr_QMixDiffCoef = xyr_TempDiffCoef
xyr_VelDiffCoef (:,:,0) = 0.0_DP
xyr_VelDiffCoef (:,:,kmax) = 0.0_DP
xyr_TempDiffCoef(:,:,0) = 0.0_DP
xyr_TempDiffCoef(:,:,kmax) = 0.0_DP
xyr_QMixDiffCoef(:,:,0) = 0.0_DP
xyr_QMixDiffCoef(:,:,kmax) = 0.0_DP
! Calculation of turbulent kinetic energy
! Turbulent kinetic energy is diagnosed.
xyr_TurKinEne = MYConstB1 * xyr_MixLength**2 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm * xyr_DVelDz**2 / 2.0_DP
xyr_TurKinEne(:,:,0 ) = 0.0_DP
xyr_TurKinEne(:,:,kmax) = 0.0_DP
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'TurKinEne', xyr_TurKinEne )
end subroutine VDiffCoefficient
| Subroutine : | |||
| 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)
| ||
| xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | 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)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated.
subroutine VDiffusionCalcFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! $ u $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! $ v $ . 南北風速. Northward wind
real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q $ . 質量混合比. Mass mixing ratio
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
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):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 質量フラックス.
! Mass flux of compositions
! 作業変数
! Work variables
!
real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP) :: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass (composition)
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 輸送係数の計算
! Calculate transfer coefficient
!
xyr_VelTransCoef (:,:,0) = 0.0_DP
xyr_VelTransCoef (:,:,kmax) = 0.0_DP
xyr_TempTransCoef(:,:,0) = 0.0_DP
xyr_TempTransCoef(:,:,kmax) = 0.0_DP
xyr_QMixTransCoef(:,:,0) = 0.0_DP
xyr_QMixTransCoef(:,:,kmax) = 0.0_DP
do k = 1, kmax-1
xyr_VelTransCoef(:,:,k) = xyr_VelDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
xyr_TempTransCoef(:,:,k) = xyr_TempDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
xyr_QMixTransCoef(:,:,k) = xyr_QMixDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
end do
! フラックスの計算
! Calculate fluxes
!
xyr_MomFluxX(:,:,0) = 0.0_DP
xyr_MomFluxX(:,:,kmax) = 0.0_DP
xyr_MomFluxY(:,:,0) = 0.0_DP
xyr_MomFluxY(:,:,kmax) = 0.0_DP
xyr_HeatFlux(:,:,0) = 0.0_DP
xyr_HeatFlux(:,:,kmax) = 0.0_DP
do n = 1, ncmax
xyrf_QMixFlux(:,:,0,n) = 0.0_DP
xyrf_QMixFlux(:,:,kmax,n) = 0.0_DP
end do
do k = 1, kmax-1
xyr_MomFluxX(:,:,k) = - xyr_VelTransCoef(:,:,k) * ( xyz_U(:,:,k+1) - xyz_U(:,:,k) )
xyr_MomFluxY(:,:,k) = - xyr_VelTransCoef(:,:,k) * ( xyz_V(:,:,k+1) - xyz_V(:,:,k) )
xyr_HeatFlux(:,:,k) = - CpDry * xyr_TempTransCoef(:,:,k) * xyr_Exner(:,:,k) * ( xyz_Temp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_Temp(:,:,k) / xyz_Exner(:,:,k) )
end do
do n = 1, ncmax
do k = 1, kmax-1
xyrf_QMixFlux(:,:,k,n) = - xyr_QMixTransCoef(:,:,k) * ( xyzf_QMix(:,:,k+1,n) - xyzf_QMix(:,:,k,n) )
end do
end do
end subroutine VDiffusionCalcFlux
| Subroutine : | |||
| z_U(1:kmax) : | real(DP), intent(in)
| ||
| z_V(1:kmax) : | real(DP), intent(in)
| ||
| zf_QMix(1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| z_Temp(1:kmax) : | real(DP), intent(in)
| ||
| r_VirTemp(0:kmax) : | real(DP), intent(in)
| ||
| r_Press(0:kmax) : | real(DP), intent(in)
| ||
| z_Height(1:kmax) : | real(DP), intent(in)
| ||
| z_Exner(1:kmax) : | real(DP), intent(in)
| ||
| r_Exner(0:kmax) : | real(DP), intent(in)
| ||
| r_VelDiffCoef(0:kmax) : | real(DP), intent(in)
| ||
| r_TempDiffCoef(0:kmax) : | real(DP), intent(in)
| ||
| r_QMixDiffCoef(0:kmax) : | real(DP), intent(in)
| ||
| r_MomFluxX(0:kmax) : | real(DP), intent(out)
| ||
| r_MomFluxY(0:kmax) : | real(DP), intent(out)
| ||
| r_HeatFlux(0:kmax) : | real(DP), intent(out)
| ||
| rf_QMixFlux(0:kmax, 1:ncmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated.
subroutine VDiffusionCalcFlux1D( z_U, z_V, zf_QMix, z_Temp, r_VirTemp, r_Press, z_Height, z_Exner, r_Exner, r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: z_U (1:kmax)
! $ u $ . 東西風速. Eastward wind
real(DP), intent(in):: z_V (1:kmax)
! $ v $ . 南北風速. Northward wind
real(DP), intent(in):: zf_QMix(1:kmax, 1:ncmax)
! $ q $ . 質量混合比. Mass mixing ratio
real(DP), intent(in):: z_Temp (1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: r_VirTemp (0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: r_Press (0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: z_Height (1:kmax)
! 高度 (整数レベル).
! Height (full level)
real(DP), intent(in):: z_Exner (1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: r_Exner (0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: r_VelDiffCoef (0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: r_TempDiffCoef (0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: r_QMixDiffCoef (0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
real(DP), intent(out):: r_MomFluxX (0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(out):: r_MomFluxY (0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(out):: r_HeatFlux (0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(out):: rf_QMixFlux(0:kmax, 1:ncmax)
! 質量フラックス.
! Mass flux of compositions
! 作業変数
! Work variables
!
real(DP) :: r_VelTransCoef (0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP) :: r_TempTransCoef(0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP) :: r_QMixTransCoef(0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass (composition)
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 輸送係数の計算
! Calculate transfer coefficient
!
r_VelTransCoef (0) = 0.0_DP
r_VelTransCoef (kmax) = 0.0_DP
r_TempTransCoef(0) = 0.0_DP
r_TempTransCoef(kmax) = 0.0_DP
r_QMixTransCoef(0) = 0.0_DP
r_QMixTransCoef(kmax) = 0.0_DP
do k = 1, kmax-1
r_VelTransCoef(k) = r_VelDiffCoef(k) * r_Press(k) / ( GasRDry * r_VirTemp(k) ) / ( z_Height(k+1) - z_Height(k) )
r_TempTransCoef(k) = r_TempDiffCoef(k) * r_Press(k) / ( GasRDry * r_VirTemp(k) ) / ( z_Height(k+1) - z_Height(k) )
r_QMixTransCoef(k) = r_QMixDiffCoef(k) * r_Press(k) / ( GasRDry * r_VirTemp(k) ) / ( z_Height(k+1) - z_Height(k) )
end do
! フラックスの計算
! Calculate fluxes
!
r_MomFluxX(0) = 0.0_DP
r_MomFluxX(kmax) = 0.0_DP
r_MomFluxY(0) = 0.0_DP
r_MomFluxY(kmax) = 0.0_DP
r_HeatFlux(0) = 0.0_DP
r_HeatFlux(kmax) = 0.0_DP
do n = 1, ncmax
rf_QMixFlux(0,n) = 0.0_DP
rf_QMixFlux(kmax,n) = 0.0_DP
end do
do k = 1, kmax-1
r_MomFluxX(k) = - r_VelTransCoef(k) * ( z_U(k+1) - z_U(k) )
r_MomFluxY(k) = - r_VelTransCoef(k) * ( z_V(k+1) - z_V(k) )
r_HeatFlux(k) = - CpDry * r_TempTransCoef(k) * r_Exner(k) * ( z_Temp(k+1) / z_Exner(k+1) - z_Temp(k) / z_Exner(k) )
end do
do n = 1, ncmax
do k = 1, kmax-1
rf_QMixFlux(k,n) = - r_QMixTransCoef(k) * ( zf_QMix(k+1,n) - zf_QMix(k,n) )
end do
end do
end subroutine VDiffusionCalcFlux1D
| Subroutine : | |||
| z_U(1:kmax) : | real(DP), intent(in)
| ||
| z_V(1:kmax) : | real(DP), intent(in)
| ||
| zf_QMix(1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| z_Temp(1:kmax) : | real(DP), intent(in)
| ||
| r_Temp(0:kmax) : | real(DP), intent(in)
| ||
| z_VirTemp(1:kmax) : | real(DP), intent(in)
| ||
| r_VirTemp(0:kmax) : | real(DP), intent(in)
| ||
| r_Press(0:kmax) : | real(DP), intent(in)
| ||
| SurfHeight : | real(DP), intent(in)
| ||
| z_Height(1:kmax) : | real(DP), intent(in)
| ||
| r_Height(0:kmax) : | real(DP), intent(in)
| ||
| z_Exner(1:kmax) : | real(DP), intent(in)
| ||
| r_Exner(0:kmax) : | real(DP), intent(in)
| ||
| z_TurKinEne(1:kmax) : | real(DP), intent(in)
| ||
| SurfMomFluxX : | real(DP), intent(in)
| ||
| SurfMomFluxY : | real(DP), intent(in)
| ||
| r_MomFluxX(0:kmax) : | real(DP), intent(out)
| ||
| r_MomFluxY(0:kmax) : | real(DP), intent(out)
| ||
| r_HeatFlux(0:kmax) : | real(DP), intent(out)
| ||
| rf_QMixFlux(0:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| r_VelDiffCoef(0:kmax) : | real(DP), intent(out)
| ||
| r_TempDiffCoef(0:kmax) : | real(DP), intent(out)
| ||
| r_QMixDiffCoef(0:kmax) : | real(DP), intent(out)
| ||
| z_DTurKinEneDt(1:kmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated by use of MY2.5 model.
subroutine VDiffusionMY251D( z_U, z_V, zf_QMix, z_Temp, r_Temp, z_VirTemp, r_VirTemp, r_Press, SurfHeight, z_Height, r_Height, z_Exner, r_Exner, z_TurKinEne, SurfMomFluxX, SurfMomFluxY, r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux, r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, z_DTurKinEneDt )
!
! 鉛直拡散フラックスを計算します.
!
! Vertical diffusion flux is calculated by use of MY2.5 model.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: FKarm, Grav, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: z_U (1:kmax)
! $ u $ . 東西風速. Eastward wind
real(DP), intent(in):: z_V (1:kmax)
! $ v $ . 南北風速. Northward wind
real(DP), intent(in):: zf_QMix(1:kmax, 1:ncmax)
! $ q $ . 質量混合比. Mass mixing ratio
real(DP), intent(in):: z_Temp (1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: r_Temp (0:kmax)
! $ \hat{T} $ . 温度 (半整数レベル).
! Temperature (half level)
real(DP), intent(in):: z_VirTemp (1:kmax)
! $ T_v $ . 仮温度. Virtual temperature
real(DP), intent(in):: r_VirTemp (0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: r_Press (0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: SurfHeight
! $ z_s $ . 地表面高度.
! Surface height.
real(DP), intent(in):: z_Height (1:kmax)
! 高度 (整数レベル).
! Height (full level)
real(DP), intent(in):: r_Height (0:kmax)
! 高度 (半整数レベル).
! Height (half level)
real(DP), intent(in):: z_Exner (1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: r_Exner (0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: z_TurKinEne(1:kmax)
!
! Turbulent kinetic energy (m2 s-2)
real(DP), intent(in):: SurfMomFluxX
!
! Eastward momentum flux at surface
real(DP), intent(in):: SurfMomFluxY
!
! Northward momentum flux at surface
real(DP), intent(out):: r_MomFluxX (0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(out):: r_MomFluxY (0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(out):: r_HeatFlux (0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(out):: rf_QMixFlux(0:kmax, 1:ncmax)
! 質量フラックス.
! Mass flux of compositions
real(DP), intent(out):: r_VelDiffCoef (0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: r_TempDiffCoef(0:kmax)
! 拡散係数：温度.
! Diffusion coefficient: temperature
real(DP), intent(out):: r_QMixDiffCoef(0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
real(DP), intent(out):: z_DTurKinEneDt (1:kmax)
!
! Tendency of turbulent kinetic energy
! 作業変数
! Work variables
!
real(DP) :: z_MixLength(1:kmax)
! 混合距離.
! Mixing length
real(DP) :: z_DVelDzSq(1:kmax)
!
! Vertical shear squared (s-2)
real(DP) :: z_StatStab(1:kmax)
!
! Static stability (s-2)
real(DP) :: GhMin
!
! Minimum of G_h
real(DP) :: GhMax
!
! Maximum of G_h
real(DP) :: z_Gm(1:kmax)
!
! G_m
real(DP) :: z_Gh(1:kmax)
!
! G_h
real(DP) :: z_Sm(1:kmax)
!
! S_M
real(DP) :: z_Sh(1:kmax)
!
! S_h
real(DP), parameter :: Stke = 0.2_DP
!
! S_{TKE} = 0.2
real(DP) :: z_VelDiffCoef (1:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP) :: z_TempDiffCoef (1:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP) :: r_TurKinEneDiffCoef (0:kmax)
!
! Diffusion coefficient: turbulent kinetic energy
real(DP) :: z_TurKinEneDiffCoef (1:kmax)
!
! Diffusion coefficient: turbulent kinetic energy
real(DP) :: r_TurKinEneTransCoef(0:kmax)
!
! Transfer coefficient: turbulent kinetic energy
real(DP) :: r_TurKinEneFlux(0:kmax)
!
! Turbulent energy flux
real(DP) :: z_CShe1(1:kmax)
real(DP) :: z_CShe2(1:kmax)
real(DP) :: z_CBuo1(1:kmax)
real(DP) :: z_CBuo2(1:kmax)
real(DP) :: z_CDis1(1:kmax)
real(DP) :: z_CDis2(1:kmax)
real(DP) :: z_TurKinEneProShear(1:kmax)
real(DP) :: z_TurKinEneProBuoya(1:kmax)
real(DP) :: FricVelSq
real(DP) :: TurKinEneAtLB
real(DP) :: za_TurKinEneMtx(1:kmax, -1:1)
!
! Implicit matrix for turbulent kinetic energy
real(DP) :: z_TurKinEneVec(1:kmax)
!
! Implicit vector for turbulent kinetic energy
real(DP) :: aaza_TurKinEneLUMtx (1:1, 1:1, 1:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP) :: aaz_DelTurKinEneLUVec(1:1, 1:1, 1:kmax)
!
! Tendency of turbulent kinetic energy
real(DP) :: z_TurKinEneDiss(1:kmax)
!
! Dissipation rate of turbulent kinetic energy (m2 s-3)
real(DP) :: z_TurKinEneNonZero(1:kmax)
!
! Turbulent kinetic energy with offset (m2 s-2)
real(DP), parameter :: TurKinEneOffset = ( 1.0e-3_DP )**2 / 2.0_DP
logical :: FlagReCalc
!
! Flag for recalculation
logical :: a_FlagReCalcLocal (1)
logical :: a_FlagReCalcGlobal(1)
integer :: iloop
integer :: nloop
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. vdiffusion_my_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 turbulent kinetic energy with offset
!
z_TurKinEneNonZero = z_TurKinEne + TurKinEneOffset
!
! Calculation of vertical shear squared
do k = 1, kmax
if ( k == 1 ) then
z_DVelDzSq(k) = ( ( z_U(k+1) - z_U(k ) )**2 + ( z_V(k+1) - z_V(k ) )**2 ) / ( z_Height(k+1) - z_Height(k ) )**2
else if ( k == kmax ) then
z_DVelDzSq(k) = ( ( z_U(k ) - z_U(k-1) )**2 + ( z_V(k ) - z_V(k-1) )**2 ) / ( z_Height(k ) - z_Height(k-1) )**2
else
z_DVelDzSq(k) = ( ( z_U(k+1) - z_U(k-1) )**2 + ( z_V(k+1) - z_V(k-1) )**2 ) / ( z_Height(k+1) - z_Height(k-1) )**2
end if
end do
! Calculation of static stability
do k = 1, kmax
if ( k == 1 ) then
z_StatStab(k) = Grav / ( z_VirTemp(k) / z_Exner(k) ) * ( z_VirTemp(k+1) / z_Exner(k+1) - z_VirTemp(k ) / z_Exner(k ) ) / ( z_Height(k+1) - z_Height(k ) )
else if ( k == kmax ) then
z_StatStab(k) = Grav / ( z_VirTemp(k) / z_Exner(k) ) * ( z_VirTemp(k ) / z_Exner(k ) - z_VirTemp(k-1) / z_Exner(k-1) ) / ( z_Height(k ) - z_Height(k-1) )
else
z_StatStab(k) = Grav / ( z_VirTemp(k) / z_Exner(k) ) * ( z_VirTemp(k+1) / z_Exner(k+1) - z_VirTemp(k-1) / z_Exner(k-1) ) / ( z_Height(k+1) - z_Height(k-1) )
end if
end do
! 混合距離の算出
! Calculate mixing length
!
do k = 1, kmax
z_MixLength(k) = FKarm * ( z_Height(k) - SurfHeight ) / (1.0_DP + FKarm * ( z_Height(k) - SurfHeight ) / MixLengthMax )
end do
! Limit mixing length (Galperin et al., 1988) and avoid zero
z_MixLength = min( z_MixLength, 0.53_DP * sqrt( 2.0_DP * z_TurKinEneNonZero / max( z_StatStab, 1.0e-10_DP ) ) ) + 1.0e-10_DP
z_Gh = - z_MixLength**2 / ( 2.0_DP * z_TurKinEneNonZero ) * z_StatStab
! Actually, xyz_Gm is not used below.
z_Gm = z_MixLength**2 / ( 2.0_DP * z_TurKinEneNonZero ) * z_DVelDzSq
! Limit Gh (Galperin et al., 1988)
GhMin = - 0.53_DP**2
GhMax = 1.0_DP / ( MYConstA2 * ( 12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) )
z_Gh = max( GhMin, min( z_Gh, GhMax ) )
z_Sh = MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) / ( 1.0_DP - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * z_Gh )
z_Sm = ( MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 - 6.0_DP * MYConstA1 / MYConstB1 ) + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) * z_Gh * z_Sh ) / ( 1.0_DP - 9.0_DP * MYConstA1 * MYConstA2 * z_Gh )
! 拡散係数の計算
! Calculation of diffusion coefficient
!
z_VelDiffCoef = z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * z_Sm
z_TempDiffCoef = z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * z_Sh
!
do k = 0, kmax
if ( ( k == 0 ) .or. ( k == kmax ) ) then
r_VelDiffCoef (k) = 0.0_DP
r_TempDiffCoef(k) = 0.0_DP
else
r_VelDiffCoef (k) = ( z_VelDiffCoef (k) + z_VelDiffCoef (k+1) ) / 2.0_DP
r_TempDiffCoef(k) = ( z_TempDiffCoef(k) + z_TempDiffCoef(k+1) ) / 2.0_DP
end if
end do
!
do k = 1, kmax-1
r_VelDiffCoef(k) = max( min( r_VelDiffCoef(k), VelDiffCoefMax ), VelDiffCoefMin )
r_TempDiffCoef(k) = max( min( r_TempDiffCoef(k), TempDiffCoefMax ), TempDiffCoefMin )
end do
!
r_QMixDiffCoef = r_TempDiffCoef
! 輸送係数とフラックスの計算
! Calculate transfer coefficient and flux
!
call VDiffusionCalcFlux1D( z_U, z_V, zf_QMix, z_Temp, r_VirTemp, r_Press, z_Height, z_Exner, r_Exner, r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux )
! Calculate tendency of turbulent kinetic energy
! Set diffusion coefficient for turbulent kinetic energy
z_TurKinEneDiffCoef = z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * Stke
!
do k = 0, kmax
if ( k == 0 ) then
r_TurKinEneDiffCoef(k) = z_TurKinEneDiffCoef(1)
else if ( k == kmax ) then
r_TurKinEneDiffCoef(k) = z_TurKinEneDiffCoef(kmax)
else
r_TurKinEneDiffCoef(k) = ( z_TurKinEneDiffCoef(k) + z_TurKinEneDiffCoef(k+1) ) / 2.0_DP
end if
end do
! Calculate turbulent kinetic energy at lower boundary
!
FricVelSq = sqrt( SurfMomFluxX**2 + SurfMomFluxY**2 ) / ( r_Press(0) / ( GasRDry * r_VirTemp(0) ) )
TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * FricVelSq
TurKinEneAtLB = TurKinEneAtLB + TurKinEneOffset
! Calculate transfer coefficient and flux of turbulent kinetic energy
!
! When transfer coefficient at lower boundary is calculated,
! diffusion coefficient at mid-point of 1st layer is used.
! In addition, transfer coefficient at upper boundary is assumed
! to be zero.
k = 0
r_TurKinEneTransCoef(k) = z_TurKinEneDiffCoef(1) * r_Press(k) / ( GasRDry * r_VirTemp(k) ) / ( z_Height(1) - SurfHeight )
do k = 1, kmax-1
r_TurKinEneTransCoef(k) = r_TurKinEneDiffCoef(k) * r_Press(k) / ( GasRDry * r_VirTemp(k) ) / ( z_Height(k+1) - z_Height(k) )
end do
k = kmax
r_TurKinEneTransCoef(k) = 0.0_DP
!
do k = 1, kmax-1
r_TurKinEneFlux(k) = - r_TurKinEneTransCoef(k) * ( z_TurKinEneNonZero(k+1) - z_TurKinEneNonZero(k) )
end do
k = 0
r_TurKinEneFlux(k) = - r_TurKinEneTransCoef(k) * ( z_TurKinEneNonZero(k+1) - TurKinEneAtLB )
k = kmax
r_TurKinEneFlux(k) = 0.0_DP
z_CShe1 = sqrt( 2.0_DP ) * z_MixLength * z_DVelDzSq * z_Sm * sqrt( z_TurKinEneNonZero )
z_CShe2 = 0.0_DP
z_CBuo1 = - sqrt( 2.0_DP ) * z_MixLength * z_StatStab * z_Sh * sqrt( z_TurKinEneNonZero )
z_CBuo2 = 0.0_DP
z_CDis1 = 2.0_DP**1.5_DP / ( MYConstB1 * z_MixLength ) * z_TurKinEneNonZero**1.5_DP
z_CDis2 = 1.5_DP * 2.0_DP**1.5_DP / ( MYConstB1 * z_MixLength ) * z_TurKinEneNonZero**0.5_DP
nloop = kmax
loop_solve : do iloop = 1, nloop
!
! Construct implicit matrix from transfer coefficient of vertical
! diffusion scheme (turbulent kinetic energy)
!
k = 1
za_TurKinEneMtx(k,-1) = 0.0_DP
za_TurKinEneMtx(k, 0) = - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) + r_TurKinEneTransCoef(k-1) + r_TurKinEneTransCoef(k ) + ( r_Press(k) - r_Press(k-1) ) / Grav * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) )
za_TurKinEneMtx(k, 1) = - r_TurKinEneTransCoef(k )
!
do k = 2, kmax-1
za_TurKinEneMtx(k,-1) = - r_TurKinEneTransCoef(k-1)
za_TurKinEneMtx(k, 0) = - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) + r_TurKinEneTransCoef(k-1) + r_TurKinEneTransCoef(k ) + ( r_Press(k) - r_Press(k-1) ) / Grav * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) )
za_TurKinEneMtx(k, 1) = - r_TurKinEneTransCoef(k )
end do
!
k = kmax
za_TurKinEneMtx(k,-1) = - r_TurKinEneTransCoef(k-1)
za_TurKinEneMtx(k, 0) = - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) + r_TurKinEneTransCoef(k-1) + r_TurKinEneTransCoef(k ) + ( r_Press(k) - r_Press(k-1) ) / Grav * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) )
za_TurKinEneMtx(k, 1) = 0.0_DP
do k = 1, kmax
z_TurKinEneVec(k) = - ( r_TurKinEneFlux(k) - r_TurKinEneFlux(k-1) ) - ( r_Press(k) - r_Press(k-1) ) / Grav * ( z_CShe1(k) + z_CBuo1(k) - z_CDis1(k) )
end do
!
! Solve simultaneous linear equations by use of LU decomposition technique
!
aaza_TurKinEneLUMtx(1,1,:,:) = za_TurKinEneMtx
!
call PhyImplLUDecomp3( aaza_TurKinEneLUMtx, 1 * 1, kmax )
aaz_DelTurKinEneLUVec(1,1,:) = z_TurKinEneVec
!
call PhyImplLUSolve3( aaz_DelTurKinEneLUVec, aaza_TurKinEneLUMtx, 1, 1 * 1 , kmax )
z_DTurKinEneDt = aaz_DelTurKinEneLUVec(1,1,:) / ( 2.0_DP * DelTime )
! Calculation of dissipation rate of turbulent kinetic energy
!
! Calculate production rate of turbulent kinetic energy
! by shear and buoyancy
z_TurKinEneProShear = z_CShe1 + z_CShe2 * z_DTurKinEneDt * 2.0_DP * DelTime
z_TurKinEneProBuoya = z_CBuo1 + z_CBuo2 * z_DTurKinEneDt * 2.0_DP * DelTime
z_TurKinEneDiss = z_CDis1 + z_CDis2 * z_DTurKinEneDt * 2.0_DP * DelTime
! Check of turbulent kinetic energy dissipation rate
! If it is negative, tendency is recalculated without dissipation.
!
FlagReCalc = .false.
do k = 1, kmax
if ( z_TurKinEneDiss(k) < 0.0_DP ) then
z_CDis1(k) = 0.0_DP
z_CDis2(k) = 0.0_DP
FlagReCalc = .true.
end if
end do
! Check convergence
a_FlagReCalcLocal = FlagReCalc
!!$ call MPIWrapperChkTrue( &
!!$ & 1, a_FlagReCalcLocal, & ! (in)
!!$ & a_FlagReCalcGlobal & ! (out)
!!$ & )
a_FlagReCalcGlobal = a_FlagReCalcLocal
if ( .not. a_FlagReCalcGlobal(1) ) exit loop_solve
end do loop_solve
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine VDiffusionMY251D
| Variable : | |||
| VelDiffCoefMax : | real(DP), save
|
| Variable : | |||
| VelDiffCoefMin : | real(DP), save
|
| Constant : | |||
| module_name = ‘vdiffusion_my‘ : | character(*), parameter
|
| Variable : | |||
| vdiffusion_my_inited = .false. : | logical, save
|
| Constant : | |||
| version = ’$Name: dcpam5-20150213 $’ // ’$Id: vdiffusion_my.f90,v 1.5 2015/01/29 12:09:17 yot Exp $’ : | character(*), parameter
|