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 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 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 : | |||
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 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 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 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 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-20131008 $’ // ’$Id: vdiffusion_my.f90,v 1.1 2013-09-16 12:01:19 yot Exp $’ : | character(*), parameter
|