| 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_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixTransCoef(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_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
!
! 鉛直拡散フラックスを計算します.
!
! 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass (composition)
! 作業変数
! 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 $
real(DP) :: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP) :: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
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.0d0
xyr_VelDiffCoef (:,:,1:kmax-1) = ConstDiffCoefM
xyr_VelDiffCoef (:,:,kmax ) = 0.0d0
xyr_TempDiffCoef(:,:,0 ) = 0.0d0
xyr_TempDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH
xyr_TempDiffCoef(:,:,kmax ) = 0.0d0
xyr_QMixDiffCoef(:,:,0 ) = 0.0d0
xyr_QMixDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH
xyr_QMixDiffCoef(:,:,kmax ) = 0.0d0
else
! バルク $ R_i $ 数算出
! Calculate bulk $ R_i $
!
xyr_DVelDz(:,:,0) = 0.
xyr_DVelDz(:,:,kmax) = 0.
xyr_BulkRiNum(:,:,0) = 0.
xyr_BulkRiNum(:,:,kmax) = 0.
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 VDiffusionTransCoefAndFlux( 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, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
! 計算時間計測一時停止
! 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.0d0
ConstDiffCoefH = 0.0d0
SquareVelMin = 0.1
BulkRiNumMin = - 100.
MixLengthMax = 300.
ShMin = 0.
SmMin = 0.
VelDiffCoefMin = 0.1
TempDiffCoefMin = 0.1
VelDiffCoefMax = 10000.
TempDiffCoefMax = 10000.
! Parameters proposed by Mellor and Yamada (1982).
!
MYConstA1 = 0.92
MYConstB1 = 16.6
MYConstA2 = 0.74
MYConstB2 = 10.1
MYConstC1 = 0.08
! 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_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixTransCoef(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_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass (composition)
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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP) :: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
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.0d-3 )**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.0d-10 ) ) ) + 1.0d-10
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 )
!!$ 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 VDiffusionTransCoefAndFlux( 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, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
! 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 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**1.5
xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**0.5
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. * 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. * 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. * 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_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixTransCoef(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_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass (composition)
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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP) :: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
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.0d-3 )**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.0d-10
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.0d100
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.0d100
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 VDiffusionTransCoefAndFlux( 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, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
! 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. * 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. * 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. * 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)
| ||
| 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_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_QMixTransCoef(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, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
!
! フラックス (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
! $ L $ [J kg-1] .
! 凝結の潜熱.
! Latent heat of condensation
! 時刻管理
! 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):: 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(in):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass of constituents
! 出力のための作業変数
! 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
!
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 )
! 風速, 温度, 比湿フラックス補正
! 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.
xyr_MomFluxXCor (:,:,kmax) = 0.
xyr_MomFluxYCor (:,:,0) = 0.
xyr_MomFluxYCor (:,:,kmax) = 0.
xyr_HeatFluxCor (:,:,0) = 0.
xyr_HeatFluxCor (:,:,kmax) = 0.
do n = 1, ncmax
xyrf_QMixFluxCor(:,:,0,n) = 0.
xyrf_QMixFluxCor(:,:,kmax,n) = 0.
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. / 3. ) - ( 2. * MYConstA1 / MYConstB1 )
Gamma2 = ( MYConstB2 / MYConstB1 ) + ( 6. * MYConstA1 / MYConstB1 )
Alpha1 = 3. * MYConstA2 * Gamma1
Alpha2 = 3. * MYConstA2 * ( Gamma1 + Gamma2 )
Beta1 = MYConstA1 * MYConstB1 * ( Gamma1 - MYConstC1 )
Beta2 = MYConstA1 * ( MYConstB1 * ( Gamma1 - MYConstC1 ) + 6. * MYConstA1 + 3. * MYConstA2 )
Beta3 = MYConstA2 * MYConstB1 * Gamma1
Beta4 = MYConstA2 * ( MYConstB1 * ( Gamma1 + Gamma2 ) - 3. * 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. * Beta2 * Beta3 * xyr_BulkRiNum ) ) / ( 2. * Beta2 )
! $ \tilde{S_h} $ と $ \tilde{S_m} $ の算出
! Calculate $ \tilde{S_h} $ and $ \tilde{S_m} $
!
xyr_Sh(:,:,kmax) = 0.
xyr_Sm(:,:,kmax) = 0.
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. - 1. * 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. + FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / MixLengthMax )
end do
! 拡散係数の算出
! Calculate diffusion constants
!
xyr_VelDiffCoef = xyr_MixLength**2 * xyr_DVelDz * sqrt ( MYConstB1 * ( 1. - xyr_FluxRiNum ) * xyr_Sm ) * xyr_Sm
xyr_TempDiffCoef = xyr_MixLength ** 2 * xyr_DVelDz * sqrt ( MYConstB1 * ( 1. - 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.
xyr_VelDiffCoef (:,:,kmax) = 0.
xyr_TempDiffCoef(:,:,0) = 0.
xyr_TempDiffCoef(:,:,kmax) = 0.
xyr_QMixDiffCoef(:,:,0) = 0.
xyr_QMixDiffCoef(:,:,kmax) = 0.
! 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_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_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixTransCoef(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 VDiffusionMY25Testing( 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_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass (composition)
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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP) :: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
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_DTurKinEneDtVDiff(0:imax-1, 1:jmax, 1:kmax)
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.0d-3 )**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
real(DP) :: xyz_TurKinEneTentative(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: DelTimeSubStep
! 実行文 ; 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.0d-10 ) ) ) + 1.0d-10
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 )
!!$ 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 VDiffusionTransCoefAndFlux( 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, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
! 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
! Calculation of vertical mixing
!
! 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. * 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. * 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. * 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_DTurKinEneDtVDiff = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime )
! Calculation tendency including source and dissipation terms
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 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**1.5
xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**0.5
! 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
xyz_TurKinEneTentative = xyz_TurKinEne
DelTimeSubStep = 0.25_DP
nloop = 2.0_DP * DelTime / DelTimeSubStep + 1
DelTimeSubStep = 2.0_DP * DelTime / dble( nloop )
nloop = kmax
loop_solve : do iloop = 1, nloop
!!$ xyz_TurKinEneProShear = &
!!$ & sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq &
!!$ & * xyz_Sm &
!!$ & * sqrt( xyz_TurKinEneTentative )
!!$ xyz_TurKinEneProBuoya = &
!!$ & - sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab &
!!$ & * xyz_Sh &
!!$ & * sqrt( xyz_TurKinEneTentative )
xyz_TurKinEneDiss = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneTentative**1.5
!!$ xyz_TurKinEneTentative = &
!!$ & xyz_TurKinEneTentative &
!!$ & + ( xyz_DTurKinEneDtVDiff &
!!$ & + xyz_TurKinEneProShear &
!!$ & + xyz_TurKinEneProBuoya &
!!$ & - xyz_TurKinEneDiss ) &
!!$ & * DelTimeSubStep
xyz_TurKinEneTentative = xyz_TurKinEneTentative + ( xyz_TurKinEneProShear + xyz_TurKinEneProBuoya - xyz_TurKinEneDiss ) * DelTimeSubStep
xyz_TurKinEneTentative = max( xyz_TurKinEneTentative, 0.0_DP )
!!$ ! 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, & ! (in)
!!$ & a_FlagReCalcGlobal & ! (out)
!!$ & )
!!$ if ( .not. a_FlagReCalcGlobal(1) ) exit loop_solve
end do loop_solve
xyz_DTurKinEneDt = ( xyz_TurKinEneTentative - xyz_TurKinEne ) / ( 2.0_DP * DelTime )
!=========================================================================
! Calculate turbulent kinetic energy with offset
!
xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset
!!$ xyz_TurKinEneNonZero = xyz_TurKinEneTentative + TurKinEneOffset
! 混合距離の算出
! 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
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 )
!!$ 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 VDiffusionTransCoefAndFlux( 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, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
! 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
! Calculation of vertical mixing
!
! 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. * 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. * 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. * 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) ) - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav * xyz_DTurKinEneDt(:,:,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 )
!!$ 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 VDiffusionMY25Testing
| 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)
| ||
| xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
|
鉛直拡散フラックスを計算します.
Vertical diffusion flux is calculated.
subroutine VDiffusionTransCoefAndFlux( 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, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
!
! 鉛直拡散フラックスを計算します.
!
! 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
real(DP), intent(out):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass (composition)
! 作業変数
! 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
! 輸送係数の計算
! Calculate transfer coefficient
!
xyr_VelTransCoef (:,:,0) = 0.
xyr_VelTransCoef (:,:,kmax) = 0.
xyr_TempTransCoef(:,:,0) = 0.
xyr_TempTransCoef(:,:,kmax) = 0.
xyr_QMixTransCoef(:,:,0) = 0.
xyr_QMixTransCoef(:,:,kmax) = 0.
do k = 1, kmax-1
!!$ xyr_VelTransCoef(:,:,k) = &
!!$ & xyr_VelDiffCoef(:,:,k) &
!!$ & * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
!!$ & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
!!$
!!$ xyr_TempTransCoef(:,:,k) = &
!!$ & xyr_TempDiffCoef(:,:,k) &
!!$ & * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
!!$ & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
!!$
!!$ xyr_QMixTransCoef(:,:,k) = &
!!$ & xyr_QMixDiffCoef(:,:,k) &
!!$ & * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
!!$ & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
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.
xyr_MomFluxX(:,:,kmax) = 0.
xyr_MomFluxY(:,:,0) = 0.
xyr_MomFluxY(:,:,kmax) = 0.
xyr_HeatFlux(:,:,0) = 0.
xyr_HeatFlux(:,:,kmax) = 0.
do n = 1, ncmax
xyrf_QMixFlux(:,:,0,n) = 0.
xyrf_QMixFlux(:,:,kmax,n) = 0.
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 VDiffusionTransCoefAndFlux
| 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-20140204-4 $’ // ’$Id: vdiffusion_my.f90,v 1.2 2013-11-17 03:18:53 yot Exp $’ : | character(*), parameter
|