| Class | dynamics_hspl_vas83 |
| In: |
dynamics/dynamics_hspl_vas83.F90
|
Note that Japanese and English are described in parallel.
力学過程を演算するモジュールです. 水平離散化にスペクトル法を, 鉛直離散化には Arakawa and Suarez (1983) を用いています.
時間積分法にはリープフロッグスキームを用いています. デフォルトでは, $ Delta t $ を大きくとるために, 重力波項に セミインプリシット法を適用しています. NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで, 重力波項をエクスプリシット法によって解くことも可能です.
This is a dynamical core module. Spectral method (for horizontal) and Arakawa and Suarez (1983) method (for vertical) are used.
Leap-frog scheme is used as time integration. By default, semi-implicit scheme is applied to gravitational terms for extension of $ Delta t $ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
| Dynamics : | 力学計算 |
| ———— : | ———— |
| Dynamics : | Calculate dynamics |
NAMELIST#dynamics_hspl_vas83_nml
| Subroutine : | |||
| xyz_UB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_QVapB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_PsB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_UN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_QVapN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_PsN(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DQVapDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_UA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_VA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_TempA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_QVapA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xy_PsA(0:imax-1, 1:jmax) : | real(DP), intent(out)
|
力学過程の演算を行い, 与えられた $ t-\Delta t $ および $ t $ の 東西風速 (xyz_UB, xyz_VN), 南北風速 (xyz_VB, xyz_VN), 温度 (xyz_TempB, xyz_TempN), 比湿 (xyz_QVapB, xyz_QVapN), 地表面気圧 (xyz_PsB, xyz_PsN) から, $ t+Delta t $ の 東西風速 (xyz_UA), 南北風速 (xyz_VA), 温度 (xyz_TempA), 比湿 (xyz_QVapA), 地表面気圧 (xyz_PsA) を返します.
別の物理プロセスによる渦度, 発散, 温度, 比湿の変化を, 力学過程の変化に足して次のステップを計算する場合には, それらの変化を xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyz_DQVapDt に与えてください.
時間積分法にはリープフロッグスキームを用いています. デフォルトでは, $ Delta t $ を大きくとるために, 重力波項に セミインプリシット法を適用しています. NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで, 重力波項をエクスプリシット法によって解くことも可能です.
Calculating dynamical core, from eastward wind (xyz_UB, xyz_UN), northward wind (xyz_VB, xyz_VN), temperature (xyz_TempB, xyz_TempN), specific humidity (xyz_QVapB, xyz_QVapN), surface pressure (xyz_PsB, xyz_PsN) at $ t-\Delta t $ and $ t $, eastward wind (xyz_UA), northward wind (xyz_VA), temperature (xyz_TempA), specific humidity (xyz_QVapA), surface pressure (xyz_PsA) are returned.
In order to add tendencies of vorticity, divergence, temperature and specific humidity by another physical process to tendencies of this dynamical process and to calculate the values at next time step, give these tendencies to "xyz_DUDt", "xyz_DVDt", "xyz_DTempDt", "xyz_DQVapDt"
Leap-frog scheme is used as time integration. By default, semi-implicit scheme is applied to gravitational terms for extension of $ Delta t $ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
subroutine Dynamics( xyz_UB, xyz_VB, xyz_TempB, xyz_QVapB, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyz_QVapN, xy_PsN, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyz_DQVapDt, xyz_UA, xyz_VA, xyz_TempA, xyz_QVapA, xy_PsA )
!
! 力学過程の演算を行い, 与えられた $ t-\Delta t $ および $ t $ の
! 東西風速 (xyz_UB, xyz_VN), 南北風速 (xyz_VB, xyz_VN),
! 温度 (xyz_TempB, xyz_TempN), 比湿 (xyz_QVapB, xyz_QVapN),
! 地表面気圧 (xyz_PsB, xyz_PsN) から,
! $ t+\Delta t $ の
! 東西風速 (xyz_UA), 南北風速 (xyz_VA), 温度 (xyz_TempA),
! 比湿 (xyz_QVapA), 地表面気圧 (xyz_PsA) を返します.
!
! 別の物理プロセスによる渦度, 発散, 温度, 比湿の変化を,
! 力学過程の変化に足して次のステップを計算する場合には,
! それらの変化を xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyz_DQVapDt
! に与えてください.
!
! 時間積分法にはリープフロッグスキームを用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! Calculating dynamical core,
! from eastward wind (xyz_UB, xyz_UN),
! northward wind (xyz_VB, xyz_VN),
! temperature (xyz_TempB, xyz_TempN),
! specific humidity (xyz_QVapB, xyz_QVapN),
! surface pressure (xyz_PsB, xyz_PsN) at $ t-\Delta t $ and $ t $,
! eastward wind (xyz_UA), northward wind (xyz_VA),
! temperature (xyz_TempA),
! specific humidity (xyz_QVapA), surface pressure (xyz_PsA)
! are returned.
!
! In order to add tendencies of vorticity, divergence,
! temperature and specific humidity by another physical process to
! tendencies of this dynamical process
! and to calculate the values at next time step,
! give these tendencies to
! "xyz_DUDt", "xyz_DVDt", "xyz_DTempDt", "xyz_DQVapDt"
!
! Leap-frog scheme is used as time integration.
! By default, semi-implicit scheme is applied to gravitational terms
! for extension of $ \Delta t $ .
! Explicit scheme can be applied to gravitational terms by changing
! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
!
! モジュール引用 ; USE statements
!
use constants, only: RPlanet, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 格子点設定
! Grid points settings
!
use gridset, only: imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
#ifdef LIB_MPI
! MPI 版 SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library MPI version, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_mpi_module, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, xy_GradLon_w => xv_GradLon_w, xy_GradLat_w => xv_GradLat_w, xya_GradLon_wa => xva_GradLon_wa, xya_GradLat_wa => xva_GradLat_wa, wa_Div_xya_xya => wa_Div_xva_xva, w_xy => w_xv, xy_w => xv_w, wa_xya => wa_xva, xya_wa => xva_wa
#else
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_module, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, xy_GradLon_w, xy_GradLat_w, xya_GradLon_wa, xya_GradLat_wa, wa_Div_xya_xya, w_xy, xy_w, wa_xya, xya_wa
#endif
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: DP ! 倍精度実数型. Double precision.
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_UB (0:imax-1, 1:jmax, 1:kmax)
! $ u (t-\Delta t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VB (0:imax-1, 1:jmax, 1:kmax)
! $ v (t-\Delta t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
! $ T (t-\Delta t) $ . 温度. Temperature
real(DP), intent(in):: xyz_QVapB (0:imax-1, 1:jmax, 1:kmax)
! $ q (t-\Delta t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsB (0:imax-1, 1:jmax)
! $ p_s (t-\Delta t) $ . 地表面気圧. Surface pressure
real(DP), intent(in):: xyz_UN (0:imax-1, 1:jmax, 1:kmax)
! $ u (t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VN (0:imax-1, 1:jmax, 1:kmax)
! $ v (t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
! $ T (t) $ . 温度. Temperature
real(DP), intent(in):: xyz_QVapN (0:imax-1, 1:jmax, 1:kmax)
! $ q (t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsN (0:imax-1, 1:jmax)
! $ p_s (t) $ . 地表面気圧. Surface pressure
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):: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{q}{t} $ . 比湿変化.
! Temperature tendency
real(DP), intent(out):: xyz_UA (0:imax-1, 1:jmax, 1:kmax)
! $ u (t+\Delta t) $ . 東西風速. Eastward wind
real(DP), intent(out):: xyz_VA (0:imax-1, 1:jmax, 1:kmax)
! $ v (t+\Delta t) $ . 南北風速. Northward wind
real(DP), intent(out):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度. Temperature
real(DP), intent(out):: xyz_QVapA (0:imax-1, 1:jmax, 1:kmax)
! $ q (t+\Delta t) $ . 比湿. Specific humidity
real(DP), intent(out):: xy_PsA (0:imax-1, 1:jmax)
! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure
! 作業変数
! Work variables
!
real(DP):: xyz_VorN (0:imax-1, 1:jmax, 1:kmax)
! $ \zeta (t) $ . 渦度. Vorticity
real(DP):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
! $ D (t) $ . 発散. Divergence
real(DP):: xy_PiN (0:imax-1, 1:jmax)
! $ \pi = \ln p_s $
real(DP):: w_PiN ((nmax+1)**2)
! $ \pi $ スペクトル
real(DP):: xy_GradLonPiN (0:imax-1, 1:jmax)
! $ \frac{1}{\cos \varphi} \DP{\pi}{\lambda} $
real(DP):: xy_GradLatPiN (0:imax-1, 1:jmax)
! $ \DP{\pi}{\varphi} $
real(DP):: xyz_UAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ u_A (t) $ . 東西運動量移流項.
! Eastward advection of momentum
real(DP):: xyz_VAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ v_A (t) $ . 南北運動量移流項.
! Northward advection of momentum
real(DP):: xyz_DTempDtN (0:imax-1, 1:jmax, 1:kmax)
! $ H (t) $ . 温度時間変化項.
! Temperature tendency
real(DP):: xyz_DQVapDtN (0:imax-1, 1:jmax, 1:kmax)
! $ R (t) $ . 比湿時間変化項.
! Specific humidity tendency
real(DP):: xyz_KinEngyN (0:imax-1, 1:jmax, 1:kmax)
! $ KE (t) $ . 運動エネルギー項.
! Kinetic energy
real(DP):: xyz_TempUAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ uT (t) $ . 温度東西移流項.
! Eastward advection of temperature
real(DP):: xyz_TempVAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ vT (t) $ . 温度南北移流項.
! Northward advection of temperature
real(DP):: xyr_SigmaDotN (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流. Vertical flow
real(DP):: xy_DPiDtN (0:imax-1, 1:jmax)
! $ Z $ . 地表面気圧時間変化項.
! Surface pressure tendency
real(DP):: xyz_QVapUAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ uq (t) $ . 比湿東西移流項.
! Eastward advection of specific humidity
real(DP):: xyz_QVapVAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ vq (t) $ . 比湿南北移流項.
! Northward advection of specific humidity
real(DP):: wz_DVorDtN ((nmax+1)**2, 1:kmax)
! $ \DD{\zeta}{t} (t) $ . 渦度変化 (スペクトル).
! Vorticity tendency (spectral)
real(DP):: wz_DDivDtN ((nmax+1)**2, 1:kmax)
! $ \DD{D}{t} (t) $ . 発散変化 (スペクトル).
! Divergence tendency (spectral)
real(DP):: wz_DTempDtN ((nmax+1)**2, 1:kmax)
! $ \DD{T}{t} (t) $ . 温度変化 (スペクトル).
! Temperature tendency (spectral)
real(DP):: wz_DQVapDtN ((nmax+1)**2, 1:kmax)
! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP):: w_DPiDtN ((nmax+1)**2)
! $ \DD{p_s}{t} (t) $ . 地表面気圧変化 (スペクトル).
! Surface pressure tendency (spectral)
real(DP):: wz_VorB ((nmax+1)**2, 1:kmax)
! $ \zeta (t-\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivB ((nmax+1)**2, 1:kmax)
! $ D (t-\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempB ((nmax+1)**2, 1:kmax)
! $ T (t-\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiB ((nmax+1)**2)
! $ \pi = \ln p_s (t-\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wz_QVapB ((nmax+1)**2, 1:kmax)
! $ q (t-\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP):: wz_VorA ((nmax+1)**2, 1:kmax)
! $ \zeta (t+\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivA ((nmax+1)**2, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempA ((nmax+1)**2, 1:kmax)
! $ T (t+\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiA ((nmax+1)**2)
! $ \pi = \ln p_s (t+\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wz_QVapA ((nmax+1)**2, 1:kmax)
! $ q (t+\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP):: wz_Psi ((nmax+1)**2, 1:kmax)
! $ \psi $ . 流線関数. Streamline function
real(DP):: wz_Chi ((nmax+1)**2, 1:kmax)
! $ \chi $ . ポテンシャル. Potential
real(DP):: wz_VorDiffA ((nmax+1)**2, 1:kmax)
! $ \mathscr{D}(\zeta) $ .
! 運動量水平拡散による渦度変化 (スペクトル).
! Vorticity tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_DivDiffA ((nmax+1)**2, 1:kmax)
! $ \mathscr{D}(D) $ .
! 運動量水平拡散による発散変化 (スペクトル).
! Divergence tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_PsiDiff ((nmax+1)**2, 1:kmax)
! 運動量水平拡散による流線関数 $ \psi $ 変化
! Streamline function tendency by
! horizontal momentum diffusion
real(DP):: wz_ChiDiff ((nmax+1)**2, 1:kmax)
! 運動量水平拡散によるポテンシャル $ \chi $ 変化
! Potential tendency by
! horizontal momentum diffusion
real(DP):: xyz_UDiff (0:imax-1, 1:jmax, 1:kmax)
! 運動量水平拡散による東西風変化.
! Eastward wind tendency by
! horizontal momentum diffusion
real(DP):: xyz_VDiff (0:imax-1, 1:jmax, 1:kmax)
! 運動量水平拡散による南北風変化.
! Northward wind tendency by
! horizontal momentum diffusion
! 外部のプロセスの時間変化項のための作業変数
! Work variables for tendency of external processes
!
real(DP):: wz_DVorDtWork ((nmax+1)**2, 1:kmax)
! $ \DD{\zeta}{t} (t) $ . 渦度変化 (スペクトル).
! Vorticity tendency (spectral)
real(DP):: wz_DDivDtWork ((nmax+1)**2, 1:kmax)
! $ \DD{D}{t} (t) $ . 発散変化 (スペクトル).
! Divergence tendency (spectral)
real(DP):: wz_DTempDtWork ((nmax+1)**2, 1:kmax)
! $ \DD{T}{t} (t) $ . 温度変化 (スペクトル).
! Temperature tendency (spectral)
real(DP):: wz_DQVapDtWork ((nmax+1)**2, 1:kmax)
! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
! エクスプリシット法のための作業変数
! Work variables for explicit scheme
!
real(DP):: xyz_exWTGPi (0:imax-1, 1:jmax, 1:kmax)
! $ \underline{W} \Dvect{T} + \Dvect{G} \pi $ .
real(DP):: xyz_exWT (0:imax-1, 1:jmax, 1:kmax)
! $ \underline{W} \Dvect{T} $ .
real(DP):: xyz_exGPi (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{G} \pi $ .
real(DP):: xyz_exHDiv (0:imax-1, 1:jmax, 1:kmax)
! $ \underline{h} \Dvect{D} $ .
integer:: k, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! 初期化
! Initialization
!
if ( .not. dynamics_hspl_vas83_inited ) call DynamicsInit
! 地表面気圧の空間変化項の計算
! Calculate spatial surface pressure tendency
!
xy_PiN = log( xy_PsN )
w_PiN = w_xy( xy_PiN )
xy_GradLonPiN = xy_GradLon_w( w_PiN )
xy_GradLatPiN = xy_GradLat_w( w_PiN )
! 風速から渦度発散の計算
! Calculate vorticity and divergence from wind velocity
!
xyz_VorN = xya_wa( wa_Div_xya_xya( xyz_VN , - xyz_UN ) / RPlanet )
xyz_DivN = xya_wa( wa_Div_xya_xya( xyz_UN , xyz_VN ) / RPlanet )
! 格子点上での非線形力学項の計算
! Calculate non-linear dynamical terms on grid points
!
call NonLinearOnGrid( xyz_UN, xyz_VN, xyz_VorN, xyz_DivN, xyz_TempN, xyz_QVapN, xy_GradLonPiN, xy_GradLatPiN, xyz_UAdvN, xyz_VAdvN, xyz_DTempDtN, xyz_DQVapDtN, xyz_KinEngyN, xyz_TempUAdvN, xyz_TempVAdvN, xyr_SigmaDotN, xy_DPiDtN, xyz_QVapUAdvN, xyz_QVapVAdvN ) ! (out)
! スペクトル時間変化項の計算
! Calculate spectral tendency terms
!
! 渦度の時間変化 (スペクトル) の計算
! Calculate vorticity tendency (spectral)
!
wz_DVorDtN = wa_Div_xya_xya( xyz_VAdvN, - xyz_UAdvN ) / RPlanet
! 発散の時間変化 (スペクトル) の計算
! Calculate divergence tendency (spectral)
!
wz_DDivDtN = wa_Div_xya_xya( xyz_UAdvN, xyz_VAdvN ) / RPlanet - wa_Lapla_wa( wa_xya(xyz_KinEngyN) ) / RPlanet**2
! 温度の時間変化 (スペクトル) の計算
! Calculate temperature tendency (spectral)
!
wz_DTempDtN = - wa_Div_xya_xya( xyz_TempUAdvN, xyz_TempVAdvN ) / RPlanet + wa_xya( xyz_DTempDtN )
! 地表面気圧の時間変化 (スペクトル) の計算
! Calculate surface pressure tendency (spectral)
!
w_DPiDtN = w_xy( xy_DPiDtN )
! 比湿の時間変化 (スペクトル) の計算
! Calculate specific humidity tendency (spectral)
!
wz_DQVapDtN = - wa_Div_xya_xya( xyz_QVapUAdvN, xyz_QVapVAdvN ) / RPlanet + wa_xya( xyz_DQVapDtN )
! エクスプリシット法を用いる際の計算
! Calculate for explicit scheme
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('explicit')
! $ \underline{W} \Dvect{T} + \Dvect{G} \pi $ の格子点値の計算
! Calculate $ \underline{W} \Dvect{T} + \Dvect{G} \pi $ on grid
!
! $ \underline{W} \Dvect{T} $ の計算
! Calculate $ \underline{W} \Dvect{T} $
!
call HydroGrid( xyz_TempN, xyz_exWT ) ! (out)
! $ \Dvect{G} \pi $ の計算
! Calculate $ \Dvect{G} \pi $
!
do k = 1, kmax
xyz_exGPi(:,:,k) = CpDry * z_TempInpolKappa(k) * z_TempAvrXY(k) * xy_PiN
enddo
! $ \underline{W} \Dvect{T} + \Dvect{G} \pi $ の計算
! Calculate $ \underline{W} \Dvect{T} + \Dvect{G} \pi $
!
xyz_exWTGPi = xyz_exWT + xyz_exGPi
! $ \underline{h} \Dvect{D} $ の格子点値の計算
! Calculate $ \underline{h} \Dvect{D} $ on grid
!
xyz_exHDiv = 0.0_DP
do k = 1, kmax
do kk = 1, kmax
xyz_exHDiv (:,:,k) = xyz_exHDiv (:,:,k) + zz_siMtxH (k,kk) * xyz_DivN (:,:,kk)
enddo
enddo
! 発散, 温度 (スペクトル) の時間変化項の修正
! Modify divergence and temperature tencency (spectral)
!
! 発散の時間変化 (スペクトル) の計算
! Calculate divergence tendency (spectral)
!
wz_DDivDtN = wz_DDivDtN - wa_Lapla_wa( wa_xya( xyz_exWTGPi ) ) / RPlanet**2
do k = 1, kmax
wz_DDivDtN(:,k) = wz_DDivDtN(:,k) - w_Lapla_w( w_Phis ) / RPlanet**2
end do
! 温度の時間変化 (スペクトル) の計算
! Calculate temperature tendency (spectral)
!
wz_DTempDtN = wz_DTempDtN - wa_xya( xyz_exHDiv )
end select
! 格子点値をスペクトル値へ ( $ t-\Delta t$ )
! Exchange grid values to spectral values ( $ t-\Delta t$ )
!
wz_VorB = wa_Div_xya_xya( xyz_VB, - xyz_UB ) / RPlanet
wz_DivB = wa_Div_xya_xya( xyz_UB, xyz_VB ) / RPlanet
wz_TempB = wa_xya( xyz_TempB )
wz_QVapB = wa_xya( xyz_QVapB )
w_PiB = w_xy( log( xy_PsB ) )
! 外部プロセスによる時間変化格子データをスペクトルデータへ変換
! Convert tendency data on grid by external processes into spectral data
!
wz_DVorDtWork = wa_Div_xya_xya( xyz_DVDt, - xyz_DUDt ) / RPlanet
wz_DDivDtWork = wa_Div_xya_xya( xyz_DUDt, xyz_DVDt ) / RPlanet
wz_DTempDtWork = wa_xya( xyz_DTempDt )
wz_DQVapDtWork = wa_xya( xyz_DQVapDt )
! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
call ImplCoef
! 時間積分
! Time integration
!
call TimeIntegration( wz_DVorDtN + wz_DVorDtWork, wz_DDivDtN + wz_DDivDtWork, wz_DTempDtN + wz_DTempDtWork, wz_DQVapDtN + wz_DQVapDtWork, w_DPiDtN, wz_VorB, wz_DivB, wz_TempB, wz_QVapB, w_PiB, wz_VorA, wz_DivA, wz_TempA, wz_QVapA, w_PiA ) ! (out)
! スペクトル値を格子点値へ ( $ t+\Delta t$ )
! Exchange spectral values to grid values ( $ t+\Delta t$ )
!
wz_Psi = wa_LaplaInv_wa( wz_VorA ) * RPlanet**2
wz_Chi = wa_LaplaInv_wa( wz_DivA ) * RPlanet**2
xyz_UA = ( xya_GradLon_wa( wz_Chi ) - xya_GradLat_wa( wz_Psi ) ) / RPlanet
xyz_VA = ( xya_GradLon_wa( wz_Psi ) + xya_GradLat_wa( wz_Chi ) ) / RPlanet
xyz_TempA = xya_wa( wz_TempA )
xyz_QVapA = xya_wa( wz_QVapA )
xy_PsA = exp( xy_w( w_PiA ) )
! 拡散による補正
! Correction by diffusion
!
! 運動量水平拡散による渦度発散の時間変化
! Vorticity and divergence tendency by horizontal diffusion of momentum
!
wz_VorDiffA = - wz_VorA * wz_DiffVorDiv
wz_DivDiffA = - wz_DivA * wz_DiffVorDiv
! 運動量水平拡散による摩擦熱補正
! Frictional thermal correction by horizontal momentum diffusion
!
wz_PsiDiff = wa_LaplaInv_wa( wz_VorDiffA ) * RPlanet**2
wz_ChiDiff = wa_LaplaInv_wa( wz_DivDiffA ) * RPlanet**2
xyz_UDiff = ( xya_GradLon_wa( wz_ChiDiff ) - xya_GradLat_wa( wz_PsiDiff ) ) / RPlanet
xyz_VDiff = ( xya_GradLon_wa( wz_PsiDiff ) + xya_GradLat_wa( wz_ChiDiff ) ) / RPlanet
xyz_TempA = xyz_TempA - ( xyz_UA * xyz_UDiff + xyz_VA * xyz_VDiff ) / CpDry * 2.0_DP * DelTime
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'SigmaDot', xyr_SigmaDotN )
call HistoryAutoPut( TimeN, 'DPiDt', xy_DPiDtN )
call DiagOutput( xyz_UA, xyz_VA, xyz_TempA, xyz_QVapA, xy_PsA ) ! (in)
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine Dynamics
| Variable : | |||
| dynamics_hspl_vas83_inited = .false. : | logical, save, public
|
| Variable : | |||
| DelTimeSave : | real(DP), save
|
| 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)
| ||
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_QVap(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_Ps(0:imax-1, 1:jmax) : | real(DP), intent(in)
|
診断量の出力を行います.
Diagnostic variables are output.
subroutine DiagOutput( xyz_U, xyz_V, xyz_Temp, xyz_QVap, xy_Ps )
!
! 診断量の出力を行います.
!
! Diagnostic variables are output.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, Grav, CpDry, LatentHeat
! $ L $ [J kg-1] .
! 凝結の潜熱.
! Latent heat of condensation
! 時刻管理
! Time control
!
use timeset, only: TimeN ! ステップ $ t $ の時刻. Time of step $ t $.
#ifdef LIB_MPI
! MPI 版 SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library MPI version, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_mpi_module, only: wa_Div_xya_xya => wa_Div_xva_xva, xya_wa => xva_wa
#else
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_module, only: wa_Div_xya_xya, xya_wa
#endif
! ヒストリデータ出力
! 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):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: xyz_QVap (0:imax-1, 1:jmax, 1:kmax)
! $ q $ . 比湿. Specific humidity
real(DP), intent(in):: xy_Ps (0:imax-1, 1:jmax)
! $ p_s $ . 地表面気圧. Surface pressure
! 作業変数
! Work variables
!
real(DP):: xyz_Vor (0:imax-1, 1:jmax, 1:kmax)
! $ \zeta $ . 渦度. Vorticity
real(DP):: xyz_Div (0:imax-1, 1:jmax, 1:kmax)
! $ D $ . 発散. Divergence
real(DP):: xy_Mass (0:imax-1, 1:jmax)
! 質量.
! Mass
real(DP):: xyz_KinEngy (0:imax-1, 1:jmax, 1:kmax)
! $ KE $ . 運動エネルギー.
! Kinetic energy
real(DP):: xyz_IntEngy (0:imax-1, 1:jmax, 1:kmax)
! $ IE $ . 内部エネルギー.
! Internal energy
real(DP):: xyz_PotEngy (0:imax-1, 1:jmax, 1:kmax)
! $ PE $ . ポテンシャルエネルギー.
! Potential energy
real(DP):: xyz_LatEngy (0:imax-1, 1:jmax, 1:kmax)
! $ LE $ . 潜熱エネルギー.
! Latent heat energy
real(DP):: xyz_TotEngy (0:imax-1, 1:jmax, 1:kmax)
! $ TE $ . 全エネルギー.
! Total energy
real(DP):: xyz_Enstro (0:imax-1, 1:jmax, 1:kmax)
! エンストロフィー.
! Enstrophy
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
! 風速から渦度発散の計算
! Calculate vorticity and divergence from wind velocity
!
xyz_Vor = xya_wa( wa_Div_xya_xya( xyz_V , - xyz_U ) / RPlanet )
xyz_Div = xya_wa( wa_Div_xya_xya( xyz_U , xyz_V ) / RPlanet )
call HistoryAutoPut( TimeN, 'Vor', xyz_Vor )
call HistoryAutoPut( TimeN, 'Div', xyz_Div )
! 質量の計算
! Calculate mass
!
xy_Mass = xy_Ps / Grav
! エネルギー, エンストロフィーの計算
! Calculate energy and enstrophy
!
call HydroGrid( xyz_Temp, xyz_PotEngy ) ! (out)
do k = 1, kmax
xyz_KinEngy(:,:,k) = ( xyz_U(:,:,k) ** 2 + xyz_V(:,:,k) ** 2 ) / 2.0_DP * xy_Mass
xyz_IntEngy(:,:,k) = CpDry * xyz_Temp(:,:,k) * xy_Mass
xyz_PotEngy(:,:,k) = xyz_PotEngy(:,:,k) * xy_Mass
xyz_LatEngy(:,:,k) = LatentHeat * xyz_QVap(:,:,k) * xy_Mass
end do
xyz_TotEngy = xyz_KinEngy + xyz_IntEngy + xyz_PotEngy + xyz_LatEngy
do k = 1, kmax
xyz_Enstro(:,:,k) = xyz_Vor(:,:,k) ** 2 * xy_Mass
end do
call HistoryAutoPut( TimeN, 'Mass', xy_Mass )
call HistoryAutoPut( TimeN, 'KinEngy', xyz_KinEngy )
call HistoryAutoPut( TimeN, 'IntEngy', xyz_IntEngy )
call HistoryAutoPut( TimeN, 'PotEngy', xyz_PotEngy )
call HistoryAutoPut( TimeN, 'LatEngy', xyz_LatEngy )
call HistoryAutoPut( TimeN, 'TotEngy', xyz_TotEngy )
call HistoryAutoPut( TimeN, 'Enstro', xyz_Enstro )
end subroutine DiagOutput
| Subroutine : |
This procedure input/output NAMELIST#dynamics_hspl_vas83_nml .
subroutine DynamicsInit
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, Omega, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 座標データ設定
! Axes data settings
!
use axesset, only: z_Sigma, r_Sigma, z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
#ifdef LIB_MPI
! MPI 版 SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library MPI version, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_mpi_module, only: xy_Lat => xv_Lat, w_xy => w_xv
#else
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_module, only: xy_Lat, w_xy
#endif
use w_module, only: rn
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoAddVariable
! gtool4 データ入力
! Gtool4 data input
!
use gtool_history, only: HistoryGet
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: STDOUT, STRING ! 文字列. Strings.
! 日付および時刻の取り扱い
! Date and time handler
!
use dc_date_types, only: DC_DIFFTIME
! 日時の差を表現するデータ型.
! Data type for difference about date and time
use dc_date, only: DCDiffTimeCreate, EvalSec
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 組み込み関数 PRESENT の拡張版関数
! Extended functions of intrinsic function "PRESENT"
!
use dc_present, only: present_and_not_empty
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 宣言文 ; Declaration statements
!
implicit none
! 地形データ (地表 $ \Phi $ ) ファイル
! File of geography data (surface $ \Phi $ )
!
character(STRING):: GeoPotentialFile
! 地形データ (地表 $ \Phi $ ) ファイル
! File of geography data (surface $ \Phi $ )
character(TOKEN):: GeoPotentialVarname
! 地形データ (地表 $ \Phi $ ) の変数名
! Variable name of geography data (surface $ \Phi $ )
! 水平拡散係数の設定のための作業変数
! Work variable for coefficient of horizontal diffusion
!
integer:: VisOrder ! 超粘性の次数. Order of hyper-viscosity
real(DP):: VisCoef ! 超粘性係数. Hyper-viscosity coefficient
real(DP):: EFoldTimeValue ! 最大波数に対する e-folding time.
! 負の値を与えると, 水平拡散係数をゼロにします.
!
! E-folding time for maximum wavenumber.
! If negative value is given,
! coefficients of horizontal diffusion become zero.
character(TOKEN):: EFoldTimeUnit
! 最大波数に対する e-folding time の単位.
! Unit of e-folding time for maximum wavenumber
real(DP):: EFoldTime ! 最大波数に対する e-folding time [単位: 秒].
! E-folding time for maximum wavenumber [Unit: sec]
type(DC_DIFFTIME):: dcdiff_efold
! NonLinearOnGrid 等で使用する係数の設定のための作業変数
! Work variable for coefficients for "NonLinearOnGrid", etc.
!
real(DP):: Kappa ! $ \kappa = R/C_p $ .
! 気体定数の定圧比熱に対する比.
! Ratio of gas constant to specific heat
! 地表ジオポテンシャル設定のための作業変数
! Work variable for surface geo-potential
!
real(DP), allocatable:: xy_Phis (:,:)
! $ \Phi_s $ . 地表ジオポテンシャル.
! Surface geo-potential
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
namelist /dynamics_hspl_vas83_nml/ TimeIntegScheme, VisOrder, EFoldTimeValue, EFoldTimeUnit, GeoPotentialFile, GeoPotentialVarname
!
! デフォルト値については初期化手続 "dynamics_hspl_vas83#DynamicsInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "dynamics_hspl_vas83#DynamicsInit" for the default values.
!
! 実行文 ; Executable statement
!
if ( dynamics_hspl_vas83_inited ) return
call InitCheck
! デフォルト値の設定
! Default values settings
!
TimeIntegScheme = 'Semi-implicit'
VisOrder = 8
EFoldTimeValue = 8640.0_DP
EFoldTimeUnit = 'sec'
GeoPotentialFile = ''
GeoPotentialVarname = ''
! 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 = dynamics_hspl_vas83_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
if ( iostat_nml == 0 ) write( STDOUT, nml = dynamics_hspl_vas83_nml )
end if
! 時間積分法のチェック
! Check time integration scheme
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
case ('explicit')
case default
call MessageNotify( 'E', module_name, '"TimeIntegScheme" must be "Semi-Implicit" or "Explicit".' )
end select
! NonLinearOnGrid 等で使用する係数の設定
! Configure coefficients for "NonLinearOnGrid", etc.
!
! コリオリパラメータの計算計算
! Calculate Coriolis parameter
!
allocate( xy_Cori (0:imax-1, 1:jmax) )
xy_Cori = 2.0_DP * Omega * sin( xy_Lat )
! 静水圧の式の係数 $ \alpha $ , $ \beta $ の計算
! Calculate coefficients $ \alpha $ and $ \beta $ in hydrostatic equation.
!
allocate( z_HydroAlpha (1:kmax) )
allocate( z_HydroBeta (1:kmax) )
Kappa = GasRDry / CpDry
do k = 1, kmax
z_HydroAlpha(k) = ( r_Sigma(k-1) / z_Sigma(k) ) ** Kappa - 1.0_DP
z_HydroBeta(k) = 1.0_DP - ( r_Sigma(k) / z_Sigma(k) ) ** Kappa
enddo
! 温度鉛直補間の係数 $ \kappa $, $ a $ , $ b $ の計算
! Calculate coefficients $ \kappa $, $ a $ , $ b $
! for interpolation of temperature
!
allocate( z_TempInpolA (1:kmax) )
allocate( z_TempInpolB (1:kmax) )
allocate( z_TempInpolKappa (1:kmax) )
do k = 1, kmax
z_TempInpolKappa(k) = ( r_Sigma(k-1) * z_HydroAlpha(k) + r_Sigma(k ) * z_HydroBeta(k) ) / z_DelSigma(k)
enddo
z_TempInpolA = 0.0_DP
do k = 2, kmax
z_TempInpolA(k) = z_HydroAlpha(k) / ( 1.0_DP - (z_Sigma(k) / z_Sigma(k-1)) ** Kappa )
end do
z_TempInpolB = 0.0_DP
do k = 1, kmax - 1
z_TempInpolB(k) = z_HydroBeta(k) / ( (z_Sigma(k) / z_Sigma(k+1) ) ** Kappa - 1.0_DP )
end do
! 平均温度 (整数レベル、半整数レベル) の計算
! Calculate average temperature about level and half level.
!
allocate( z_TempAvrXY (1:kmax) )
allocate( r_TempAvrXY (0:kmax) )
z_TempAvrXY = 300.0_DP
r_TempAvrXY = 0.0_DP
do k = 1, kmax - 1
r_TempAvrXY(k) = z_TempInpolA(k+1) * z_TempAvrXY(k+1) + z_TempInpolB( k ) * z_TempAvrXY( k )
enddo
! 水平拡散係数の設定
! Configure coefficient of horizontal diffusion
!
allocate( wz_rn ((nmax+1)**2, 1:kmax) )
allocate( wz_DiffVorDiv ((nmax+1)**2, 1:kmax) )
allocate( wz_DiffTherm ((nmax+1)**2, 1:kmax) )
do k = 1, kmax
wz_rn(:,k) = rn(:,1)
enddo
! 粘性係数の計算 (最大波数の e-folding time が EFoldTime となるように)
! Calculate viscosity coefficient
!
call DCDiffTimeCreate( dcdiff_efold, EFoldTimeValue, EFoldTimeUnit ) ! (in)
EFoldTime = EvalSec( dcdiff_efold )
if ( EFoldTimeValue > 0.0_DP ) then
VisCoef = ( (nmax*(nmax+1)) / RPlanet**2 )**(-VisOrder / 2) / EFoldTime
else
VisCoef = 0.0_DP
end if
wz_DiffTherm = VisCoef * ( (-wz_rn / RPlanet**2)**(VisOrder / 2) )
wz_DiffVorDiv = wz_DiffTherm + VisCoef * ( - (2.0_DP / RPlanet**2)**(VisOrder / 2))
! 地表ジオポテンシャルの設定
! Configure surface geo-potential
!
allocate( xy_Phis (0:imax-1, 1:jmax) )
allocate( w_Phis ((nmax+1)**2) )
if ( present_and_not_empty( GeoPotentialFile ) .and. present_and_not_empty( GeoPotentialVarname ) ) then
call HistoryGet( GeoPotentialFile, GeoPotentialVarname, xy_Phis ) ! (out)
w_Phis = w_xy( xy_Phis )
else
w_Phis = 0.0_DP
end if
! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
call ImplCoef
! ヒストリデータ出力のためのへの変数登録
! Register of variables for history data output
!
call HistoryAutoAddVariable( 'SigmaDot', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'sigma-vertical velocity', '1 s-1' )
call HistoryAutoAddVariable( 'DPiDt', (/ 'lon ', 'lat ', 'time' /), 'Pi (log Ps) tendency)', 'Pa s-1' )
call HistoryAutoAddVariable( 'Vor', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'vorticity', 's-1' )
call HistoryAutoAddVariable( 'Div', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'divergence', 's-1' )
call HistoryAutoAddVariable( 'Mass', (/ 'lon ', 'lat ', 'time' /), 'mass', 'kg' )
call HistoryAutoAddVariable( 'KinEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'kinetic energy', 'J' )
call HistoryAutoAddVariable( 'IntEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'internal energy', 'J' )
call HistoryAutoAddVariable( 'PotEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'potential energy', 'J' )
call HistoryAutoAddVariable( 'LatEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'latent energy', 'J' )
call HistoryAutoAddVariable( 'TotEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'total energy', 'J' )
call HistoryAutoAddVariable( 'Enstro', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'enstrophy', 'kg' )
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' TimeIntegScheme = %c', c1 = trim( TimeIntegScheme ) )
call MessageNotify( 'M', module_name, ' EFoldTime = %f [%c]', d = (/ EFoldTimeValue /), c1 = trim(EFoldTimeUnit) )
call MessageNotify( 'M', module_name, ' VisOrder = %d', i = (/ VisOrder /) )
call MessageNotify( 'M', module_name, ' VisCoef = %f', d = (/ VisCoef /) )
call MessageNotify( 'M', module_name, ' GeoPotentialFile = %c', c1 = trim( GeoPotentialFile ) )
call MessageNotify( 'M', module_name, ' GeoPotentialVarname = %c', c1 = trim( GeoPotentialVarname ) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
dynamics_hspl_vas83_inited = .true.
end subroutine DynamicsInit
| Subroutine : | |||
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_Phi(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
格子点データである温度 $ T $ から, 静水圧の式を用いて 格子点データのジオポテンシャル高度 $ Phi $ を求めます.
subroutine HydroGrid( xyz_Temp, xyz_Phi )
!
! 格子点データである温度 $ T $ から, 静水圧の式を用いて
! 格子点データのジオポテンシャル高度 $ \Phi $ を求めます.
!
! モジュール引用 ; USE statements
!
use constants, only: CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(out):: xyz_Phi (0:imax-1, 1:jmax, 1:kmax)
! $ \Phi $ . ジオポテンシャル高度.
! Getpotential height
! 作業変数
! Work variables
!
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
xyz_Phi(:,:,1) = CpDry * z_HydroAlpha(1) * xyz_Temp(:,:,1)
do k = 2, kmax
xyz_Phi(:,:,k) = xyz_Phi(:,:,k-1) + CpDry * z_HydroAlpha(k) * xyz_Temp(:,:,k) + CpDry * z_HydroBeta(k-1) * xyz_Temp(:,:,k-1)
enddo
end subroutine HydroGrid
| Subroutine : |
TimeIntegration で使用する係数の設定.
Configure coefficients for "TimeIntegration"
subroutine ImplCoef
!
! TimeIntegration で使用する係数の設定.
!
! Configure coefficients for "TimeIntegration"
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 座標データ設定
! Axes data settings
!
use axesset, only: r_Sigma, z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_module, only: l_nm
! LU 分解法により連立 1 次方程式を解くための関数
! Functions to solve linear simultaneous equation by LU decomposition
!
use lumatrix, only: LUDecomp
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! デバッグ用ユーティリティ
! Utilities for debug
!
use dc_trace, only: DbgMessage, BeginSub, EndSub, Debug
! 宣言文 ; Declaration statements
!
implicit none
! TimeIntegration 等で使用する係数の設定のための作業変数
! Work variable for coefficients for "TimeIntegration", etc.
!
real(DP), allocatable:: zz_siMtxW (:,:)
! $ W $ .
! 発散の式での線形重力波項の効果による温度補正の係数.
! Coefficient for correction of temperature
! by effort of linear gravitational terms
real(DP), allocatable:: zz_siMtxQ (:,:)
! $ Q = \DD{T}{\sigma} $
real(DP), allocatable:: zz_siMtxS (:,:)
! $ S = \DD{\sigma}{t} $
real(DP), allocatable:: zz_siMtxQS (:,:)
! $ QS $ .
! この変数は $ \sigma $ 移流の効果に相当.
! This variable corresponds to effort of $ \sigma $ advection
real(DP), allocatable:: zz_siMtxR (:,:)
! $ R $ .
! 発散と掛け合わせることで気圧変化の効果となる.
! Product R and divergence become effort of
! surface pressure tendency.
! $ RD = \kappa T
! (\DD{\pi}{t} + \Dinv{\sigma}\DD{\sigma}{t}) $ .
integer, allocatable:: nmo (:,:,:)
! スペクトルの添字順番.
! Spectral subscript expression
real(DP), allocatable:: zz_siMtxM (:,:)
! 行列 $ \underline{M} $.
! Matrix $ \underline{M} $
integer, allocatable:: z_siMtxPivWork(:)
! 行列のピボット作成の作業変数.
! Work variable for pivot
real(DP):: flapla ! $ \nabla_{\sigma}^{2} $
integer:: k, l, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: mmax ! 最大東西波数.
! Maximum truncated eastward wavenumber
integer:: lmax ! 最大南北波数.
! Maximum truncated northward wavenumber
integer:: n, m, nm, mxnm
! 波数方向に回る DO ループ用作業変数
! Work variables for DO loop in wavenumber direction
logical:: lmax_err ! 最大南北波数に関するエラーフラグ.
! Error flag for maximum truncated northward wavenumber
! 実行文 ; Executable statement
!
! $ \Delta t $ [s] のチェック・保存
! Check and save $ \Delta t $ [s]
!
if ( .not. dynamics_hspl_vas83_inited ) then
DelTimeSave = DelTime
else
if ( DelTimeSave == DelTime ) then
return
else
DelTimeSave = DelTime
end if
end if
call DbgMessage( '%c:%c: (DelTime=%f) coefficients for "TimeIntegration" is generated. ', c1 = module_name, c2 = 'ImplCoef', d = (/ DelTime /) )
! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
if ( .not. allocated( z_siMtxG ) ) allocate( z_siMtxG (1:kmax) )
if ( .not. allocated( zz_siMtxH ) ) allocate( zz_siMtxH (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxWH ) ) allocate( zz_siMtxWH (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxGCt) ) allocate( zz_siMtxGCt (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxW ) ) allocate( zz_siMtxW (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxQ ) ) allocate( zz_siMtxQ (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxS ) ) allocate( zz_siMtxS (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxQS ) ) allocate( zz_siMtxQS (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxR ) ) allocate( zz_siMtxR (1:kmax, 1:kmax) )
z_siMtxG = CpDry * z_TempInpolKappa * z_TempAvrXY
do k = 1, kmax
do l = 1, kmax
zz_siMtxGCt(k,l) = z_siMtxG(k) * z_DelSigma(l)
end do
end do
zz_siMtxW = 0.0_DP
do k = 1, kmax
do l = 1, k
zz_siMtxW(k,l) = CpDry * z_HydroAlpha(l)
enddo
do l = 1, k-1
zz_siMtxW(k,l) = zz_siMtxW(k,l) + CpDry * z_HydroBeta(l)
enddo
enddo
zz_siMtxS = 0.0_DP
do k = 1, kmax
do l = 1, kmax
zz_siMtxS(k,l) = r_Sigma(k-1) * z_DelSigma(l)
enddo
do l = k, kmax
zz_siMtxS(k,l) = zz_siMtxS(k,l) - z_DelSigma(l)
enddo
enddo
zz_siMtxQ = 0.0_DP
do k = 1, kmax
zz_siMtxQ(k,k) = ( r_TempAvrXY(k-1) - z_TempAvrXY(k) ) / z_DelSigma(k)
enddo
do k = 1, kmax-1
zz_siMtxQ(k,k+1) = ( z_TempAvrXY(k) - r_TempAvrXY(k) ) / z_DelSigma(k)
enddo
zz_siMtxQS = 0.0_DP
zz_siMtxQS = matmul(zz_siMtxQ, zz_siMtxS)
zz_siMtxR = 0.0_DP
do k = 1, kmax
do l = k, kmax
zz_siMtxR(k,l) = - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(l) * z_TempAvrXY(k)
enddo
do l = k + 1, kmax
zz_siMtxR(k,l) = zz_siMtxR(k,l) - z_HydroBeta(k) / z_DelSigma(k) * z_DelSigma(l) * z_TempAvrXY(k)
enddo
enddo
zz_siMtxH = 0.0_DP
zz_siMtxH = zz_siMtxQS - zz_siMtxR
zz_siMtxWH = 0.0_DP
zz_siMtxWH = matmul(zz_siMtxW, zz_siMtxH)
if ( .not. allocated(nmo ) ) allocate( nmo (1:2, 0:nmax, 0:nmax) )
if ( .not. allocated(zz_siMtxM ) ) allocate( zz_siMtxM (1:kmax, 1:kmax) )
if ( .not. allocated(z_siMtxPivWork) ) allocate( z_siMtxPivWork(1:kmax) )
if ( .not. allocated(wzz_siMtxLU ) ) allocate( wzz_siMtxLU ((nmax+1)**2, 1:kmax, 1:kmax) )
if ( .not. allocated(wz_siMtxPiv ) ) allocate( wz_siMtxPiv ((nmax+1)**2, 1:kmax) )
mmax = nmax
lmax = nmax
mxnm = 0
! スペクトル添字順序の取り出し
! Fetch spectral subscript expression
!
nmo = 0
do l = 0, lmax
do m = 0, min(mmax, nmax-l)
nmo(1,m,l) = l_nm(m+l, m)
nmo(2,m,l) = l_nm(m+l, -m)
enddo
enddo
Loop_n: do n = 0, nmax
flapla = - real(n) * real(n+1) / RPlanet**2
! スペクトル添字順序の取り出し
! Fetch spectral subscript expression
!
lmax_err = .true.
do m = 0, min(n,mmax)
if ( n-m <= lmax ) then
nm = nmo(1,m,n-m)
lmax_err = .false.
exit
endif
end do
if ( lmax_err ) then
call MessageNotify( 'E', module_name, 'n-m=<%d> must be less than or equal to lmax=<%d>.', i = (/ n-m, lmax /) )
end if
! 行列 $ \underline{M} $ の計算
! Calculate matrix $ \underline{M} $
!
do k = 1, kmax
do kk = 1, kmax
zz_siMtxM ( k,kk ) = - DelTime**2 * flapla * ( zz_siMtxWH( k,kk ) + zz_siMtxGCt( k,kk ) * ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(nm,1) ) )
if ( k == kk ) then
zz_siMtxM ( k,kk ) = zz_siMtxM ( k,kk ) + ( 1.0_DP + DelTime * 2.0_DP * wz_DiffVorDiv(nm,1) ) * ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(nm,1) )
endif
end do
end do
! LU 行列計算
! LU matrix calculation
!
call LUDecomp( zz_siMtxM, z_siMtxPivWork ) ! (out)
! ダミー値の代入. (位置 kmax はまだ未定義なため).
! Dummy value is subtituted. (Because position kmax is undefined yet).
!
z_siMtxPivWork(kmax) = 0
!!$ write(*,*) 'n= ', n
!!$
!!$ do k = 1, kmax
!!$ do kk = 1, kmax
!!$ write(*,*) 'zz_siMtxM(', k+1, ',', kk+1, ')= ', zz_siMtxM(k,kk)
!!$ end do
!!$ end do
!!$ do k = 1, kmax
!!$ write(*,*) 'z_siMtxPivWork(', k, ')= ', z_siMtxPivWork(k)
!!$ end do
! 配列の詰め替え
! Repack matrices
!
do m = 0, mmax
l = n - m
if ( ( l >= 0 ) .and. ( l <= lmax ) ) then
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxLU ( nmo(1,m,l),k,kk ) = zz_siMtxM ( k,kk )
wzz_siMtxLU ( nmo(2,m,l),k,kk ) = zz_siMtxM ( k,kk )
end do
wz_siMtxPiv ( nmo(1,m,l),k ) = z_siMtxPivWork ( k )
wz_siMtxPiv ( nmo(2,m,l),k ) = z_siMtxPivWork ( k )
mxnm = max( mxnm, nmo(1,m,l) )
mxnm = max( mxnm, nmo(2,m,l) )
end do
endif
end do
end do Loop_n
do nm = mxnm+1, (nmax+1)**2
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxLU ( nm,k,kk ) = zz_siMtxM ( k,kk )
end do
wz_siMtxPiv ( nm,k ) = z_siMtxPivWork ( k )
end do
end do
end subroutine ImplCoef
| Subroutine : |
依存モジュールの初期化チェック
Check initialization of dependency modules
subroutine InitCheck
!
! 依存モジュールの初期化チェック
!
! Check initialization of dependency modules
! モジュール引用 ; USE statements
!
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_util_inited
! 格子点設定
! Grid points settings
!
use gridset, only: gridset_inited
! 物理定数設定
! Physical constants settings
!
use constants, only: constants_inited
! 座標データ設定
! Axes data settings
!
use axesset, only: axesset_inited
! 時刻管理
! Time control
!
use timeset, only: timeset_inited
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 実行文 ; Executable statement
!
if ( .not. namelist_util_inited ) call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )
if ( .not. gridset_inited ) call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )
if ( .not. constants_inited ) call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )
if ( .not. axesset_inited ) call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )
if ( .not. timeset_inited ) call MessageNotify( 'E', module_name, '"timeset" module is not initialized.' )
end subroutine InitCheck
| 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)
| ||
| xyz_Vor(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_Div(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_QVap(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_GradLonPi(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_GradLatPi(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_UAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_VAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_DQVapDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_KinEngy(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_TempUAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_TempVAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyr_SigmaDot(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xy_DPiDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xyz_QVapUAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_QVapVAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
非線形項 (非重力波項) を格子点上で計算します.
Non-linear terms (non gravitational terms) are calculated on grid points
subroutine NonLinearOnGrid( xyz_U, xyz_V, xyz_Vor, xyz_Div, xyz_Temp, xyz_QVap, xy_GradLonPi, xy_GradLatPi, xyz_UAdv, xyz_VAdv, xyz_DTempDt, xyz_DQVapDt, xyz_KinEngy, xyz_TempUAdv, xyz_TempVAdv, xyr_SigmaDot, xy_DPiDt, xyz_QVapUAdv, xyz_QVapVAdv )
!
! 非線形項 (非重力波項) を格子点上で計算します.
!
! Non-linear terms (non gravitational terms) are calculated on
! grid points
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: r_Sigma, z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, CpDry, EpsV
! $ \epsilon_v $ .
! 水蒸気分子量比.
! Molecular weight of water vapor
! 文字列操作
! Character handling
!
use dc_string, only: LChar
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):: xyz_Vor (0:imax-1, 1:jmax, 1:kmax)
! $ \zeta $ . 渦度. Vorticity
real(DP), intent(in):: xyz_Div (0:imax-1, 1:jmax, 1:kmax)
! $ D $ . 発散. Divergence
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(in):: xyz_QVap (0:imax-1, 1:jmax, 1:kmax)
! $ q $ . 比湿. Specific humidity
real(DP), intent(in):: xy_GradLonPi (0:imax-1, 1:jmax)
! $ \frac{1}{\cos \varphi} \DP{\pi}{\lambda} $
real(DP), intent(in):: xy_GradLatPi (0:imax-1, 1:jmax)
! $ \DP{\pi}{\varphi} $
real(DP), intent(out):: xyz_UAdv (0:imax-1, 1:jmax, 1:kmax)
! $ u_A $ . 東西運動量移流項.
! Eastward advection of momentum
real(DP), intent(out):: xyz_VAdv (0:imax-1, 1:jmax, 1:kmax)
! $ v_A $ . 南北運動量移流項.
! Northward advection of momentum
real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ H $ . 温度時間変化項.
! Temperature tendency
real(DP), intent(out):: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax)
! $ R $ . 比湿時間変化項.
! Specific humidity tendency
real(DP), intent(out):: xyz_KinEngy (0:imax-1, 1:jmax, 1:kmax)
! $ KE $ . 運動エネルギー項.
! Kinetic energy
real(DP), intent(out):: xyz_TempUAdv (0:imax-1, 1:jmax, 1:kmax)
! $ uT $ . 温度東西移流項.
! Eastward advection of temperature
real(DP), intent(out):: xyz_TempVAdv (0:imax-1, 1:jmax, 1:kmax)
! $ vT $ . 温度南北移流項.
! Northward advection of temperature
real(DP), intent(out):: xyr_SigmaDot (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流. Vertical flow
real(DP), intent(out):: xy_DPiDt (0:imax-1, 1:jmax)
! $ Z $ . 地表面気圧時間変化項.
! Surface pressure tendency
real(DP), intent(out):: xyz_QVapUAdv (0:imax-1, 1:jmax, 1:kmax)
! $ uq $ . 比湿東西移流項.
! Eastward advection of specific humidity
real(DP), intent(out):: xyz_QVapVAdv (0:imax-1, 1:jmax, 1:kmax)
! $ vq $ . 比湿南北移流項.
! Northward advection of specific humidity
!-----------------------------------
! 作業変数
! Work variables
real(DP):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{v} \cdot \nabla \pi $ .
! $ \pi $ の移流. Advection of $ \pi $
real(DP):: xyz_PiAdvSum (0:imax-1, 1:jmax, 1:kmax)
! $ \sum_k^K(\Dvect{v}\cdot\nabla\pi)\Delta\sigma $ .
! $ \pi $ 移流の積下げ. Integral downward of advection of $ \pi $
real(DP):: xyz_DivSum (0:imax-1, 1:jmax, 1:kmax)
! $ \sum_k^K D\Delta\sigma $ .
! 発散の積下げ. Integral downward of divergence
real(DP):: xyr_SigmaDotNonGrav (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流 (非重力波). Vertical flow (non gravitational)
real(DP):: xyz_TempEdd (0:imax-1, 1:jmax, 1:kmax)
! $ T' = T - \bar{T} $ .
! 温度の擾乱 (整数レベル). Temperature eddy (full level)
real(DP):: xyr_TempEdd (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}' $ .
! 温度の擾乱 (半整数レベル). Temperature eddy (half level)
real(DP):: xyz_TempVir (0:imax-1, 1:jmax, 1:kmax)
! $ T_v $ .
! 仮温度. Virtual temperature
real(DP):: xyz_TempVirEdd (0:imax-1, 1:jmax, 1:kmax)
! $ {T_v}' = T_v - \bar{T} $ .
! 仮温度の擾乱. Virtual temperature eddy
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
! $ \pi $ の移流, $ \pi $ 移流の積下げ, 発散の積下げの計算
! Calculate advection of $ \pi $, integral downward of advection
! of $ \pi $, integral downward of divergence
!
do k = 1, kmax
xyz_PiAdv(:,:,k) = ( xyz_U(:,:,k) * xy_GradLonPi + xyz_V(:,:,k) * xy_GradLatPi ) / RPlanet
enddo
xyz_PiAdvSum(:,:,kmax) = xyz_PiAdv(:,:,kmax) * z_DelSigma(kmax)
do k = kmax-1, 1, -1
xyz_PiAdvSum(:,:,k) = xyz_PiAdvSum(:,:,k+1) + xyz_PiAdv(:,:,k) * z_DelSigma(k)
enddo
xyz_DivSum(:,:,kmax) = xyz_Div(:,:,kmax) * z_DelSigma(kmax)
do k = kmax-1, 1, -1
xyz_DivSum(:,:,k) = xyz_DivSum(:,:,k+1) + xyz_Div(:,:,k) * z_DelSigma(k)
enddo
! 地表面気圧時間変化項 $ Z $ の計算
! Calculate surface pressure tendency $ Z $
!
xy_DPiDt = - xyz_PiAdvSum(:,:,1)
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('explicit')
xy_DPiDt = xy_DPiDt - xyz_DivSum(:,:,1)
end select
! $ \dot{\sigma} $ の計算
! Calculate $ \dot{\sigma} $
!
do k = 1, kmax-1
xyr_SigmaDot(:,:,k) = r_Sigma(k) * ( xyz_PiAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) - ( xyz_PiAdvSum(:,:,k+1) + xyz_DivSum(:,:,k+1) )
xyr_SigmaDotNonGrav(:,:,k) = r_Sigma(k) * xyz_PiAdvSum(:,:,1) - xyz_PiAdvSum(:,:,k+1)
enddo
! $ \dot{\sigma} $ の上下境界値
! $ \dot{\sigma} $ on upper and lower boundary
!
xyr_SigmaDot(:,:,0) = 0.0_DP
xyr_SigmaDot(:,:,kmax) = 0.0_DP
xyr_SigmaDotNonGrav(:,:,0) = 0.0_DP
xyr_SigmaDotNonGrav(:,:,kmax) = 0.0_DP
! 温度の擾乱 (整数レベル), 仮温度, 仮温度の擾乱の計算
! Calculate temperature eddy (full level), virtual temperature,
! virtual temperature eddy
!
do k = 1, kmax
xyz_TempVir(:,:,k) = xyz_Temp(:,:,k) * ( 1.0_DP + ( ( ( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QVap(:,:,k) ) )
xyz_TempEdd(:,:,k) = xyz_Temp(:,:,k) - z_TempAvrXY(k)
xyz_TempVirEdd(:,:,k) = xyz_TempVir(:,:,k) - z_TempAvrXY(k)
enddo
! 温度の擾乱 (半整数レベル) の計算
! Calculate temperature eddy (half level)
!
xyr_TempEdd(:,:,0) = 0.0_DP
xyr_TempEdd(:,:,kmax) = 0.0_DP
do k = 1, kmax-1
xyr_TempEdd(:,:,k) = z_TempInpolA(k+1) * xyz_Temp(:,:,k+1) + z_TempInpolB( k ) * xyz_Temp(:,:, k ) - r_TempAvrXY(k)
enddo
! 東西運動量移流項, 南北運動量移流項の計算
! Calculate advection of eastward momentum and northward momentum
!
do k = 1, kmax
xyz_UAdv(:,:,k) = ( xyz_Vor(:,:,k) + xy_Cori ) * xyz_V(:,:,k) - CpDry * z_TempInpolKappa(k) * xyz_TempVirEdd(:,:,k) * xy_GradLonPi / RPlanet
xyz_VAdv(:,:,k) = - ( xyz_Vor(:,:,k) + xy_Cori ) * xyz_U(:,:,k) - CpDry * z_TempInpolKappa(k) * xyz_TempVirEdd(:,:,k) * xy_GradLatPi / RPlanet
end do
do k = 2, kmax
xyz_UAdv(:,:,k) = xyz_UAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k-1) * ( xyz_U(:,:,k-1) - xyz_U(:,:,k) )
xyz_VAdv(:,:,k) = xyz_VAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k-1) * ( xyz_V(:,:,k-1) - xyz_V(:,:,k) )
end do
do k = 1, kmax-1
xyz_UAdv(:,:,k) = xyz_UAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k) * ( xyz_U(:,:,k) - xyz_U(:,:,k+1) )
xyz_VAdv(:,:,k) = xyz_VAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k) * ( xyz_V(:,:,k) - xyz_V(:,:,k+1) )
end do
! 運動エネルギー項 (仮温度補正含む) の計算
! Calculate kinematic energy term
! (including virtual temperature correction)
!
call HydroGrid( xyz_TempVir - xyz_Temp, xyz_KinEngy ) ! (out)
xyz_KinEngy = xyz_KinEngy + ( xyz_U**2 + xyz_V**2 ) / 2.0_DP
! 温度東西移流項, 温度南北移流項の計算
! Calculate eastward and northward advection of temperature
!
do k = 1, kmax
xyz_TempUAdv(:,:,k) = xyz_U(:,:,k) * xyz_TempEdd(:,:,k)
xyz_TempVAdv(:,:,k) = xyz_V(:,:,k) * xyz_TempEdd(:,:,k)
end do
! 温度の時間変化項 $ H $ の計算
! Calculate temperature tendency term $ H $
!
do k = 1, kmax
xyz_DTempDt(:,:,k) = xyz_TempEdd(:,:,k) * xyz_Div(:,:,k) + z_TempInpolKappa(k) * xyz_TempVir(:,:,k) * xyz_PiAdv(:,:,k) - z_HydroAlpha(k) / z_DelSigma(k) * ( xyz_TempVir(:,:,k) * xyz_PiAdvSum(:,:,k) + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) )
enddo
do k = 2, kmax
xyz_DTempDt(:,:,k) = xyz_DTempDt(:,:,k) - xyr_SigmaDot(:,:,k-1) / z_DelSigma(k) * ( xyr_TempEdd(:,:,k-1) - xyz_TempEdd(:,:,k) ) - xyr_SigmaDotNonGrav(:,:,k-1) / z_DelSigma(k) * ( r_TempAvrXY(k-1) - z_TempAvrXY(k) )
enddo
do k = 1, kmax-1
xyz_DTempDt(:,:,k) = xyz_DTempDt(:,:,k) - xyr_SigmaDot(:,:,k) / z_DelSigma(k) * ( xyz_TempEdd(:,:,k) - xyr_TempEdd(:,:,k) ) - xyr_SigmaDotNonGrav(:,:,k) / z_DelSigma(k) * ( z_TempAvrXY(k) - r_TempAvrXY(k) ) - z_HydroBeta(k) / z_DelSigma(k) * ( xyz_TempVir(:,:,k) * xyz_PiAdvSum(:,:,k+1) + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1) )
enddo
! 比湿東西移流項, 比湿南北移流項の計算
! Calculate eastward and northward advection of specific humidity
!
do k = 1, kmax
xyz_QVapUAdv(:,:,k) = xyz_U(:,:,k) * xyz_QVap(:,:,k)
xyz_QVapVAdv(:,:,k) = xyz_V(:,:,k) * xyz_QVap(:,:,k)
end do
! 比湿時間変化項 $ R $ の計算
! Calculate specific humidity tendency $ R $
!
xyz_DQVapDt = xyz_QVap * xyz_Div
do k = 2, kmax
xyz_DQVapDt(:,:,k) = xyz_DQVapDt(:,:,k) - xyr_SigmaDot(:,:,k-1) / ( 2.0_DP * z_DelSigma(k) ) * ( xyz_QVap(:,:,k-1) - xyz_QVap(:,:,k) )
enddo
do k = 1, kmax-1
xyz_DQVapDt(:,:,k) = xyz_DQVapDt(:,:,k) - xyr_SigmaDot(:,:,k) / ( 2.0_DP * z_DelSigma(k) ) * ( xyz_QVap(:,:,k) - xyz_QVap(:,:,k+1) )
enddo
end subroutine NonLinearOnGrid
| Variable : | |||
| TimeIntegScheme : | character(TOKEN), save
|
| Subroutine : | |||
| wz_DVorDtN((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| wz_DDivDtN((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| wz_DTempDtN((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| wz_DQVapDtN((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| w_DPiDtN((nmax+1)**2) : | real(DP), intent(in)
| ||
| wz_VorB((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| wz_DivB((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| wz_TempB((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| wz_QVapB((nmax+1)**2, 1:kmax) : | real(DP), intent(in)
| ||
| w_PiB((nmax+1)**2) : | real(DP), intent(in)
| ||
| wz_VorA((nmax+1)**2, 1:kmax) : | real(DP), intent(out)
| ||
| wz_DivA((nmax+1)**2, 1:kmax) : | real(DP), intent(out)
| ||
| wz_TempA((nmax+1)**2, 1:kmax) : | real(DP), intent(out)
| ||
| wz_QVapA((nmax+1)**2, 1:kmax) : | real(DP), intent(out)
| ||
| w_PiA((nmax+1)**2) : | real(DP), intent(out)
|
時間積分を行い, 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から 時刻 $ t+Delta t $ の物理量を計算します.
時間積分法にはリープフロッグスキームを用いています. ただし拡散項には前方差分を用いています. デフォルトでは, $ Delta t $ を大きくとるために, 重力波項に セミインプリシット法を適用しています. NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで, 重力波項をエクスプリシット法によって解くことも可能です.
With time integration, calculate physical values at $ t+Delta t $ from tendency at $ t $ and physical values at $ t-\Delta t $ .
Leap-frog scheme is used as time integration scheme. And forward difference is used to diffusion terms. By default, semi-implicit scheme is applied to gravitational terms for extension of $ Delta t $ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
subroutine TimeIntegration( wz_DVorDtN, wz_DDivDtN, wz_DTempDtN, wz_DQVapDtN, w_DPiDtN, wz_VorB, wz_DivB, wz_TempB, wz_QVapB, w_PiB, wz_VorA, wz_DivA, wz_TempA, wz_QVapA, w_PiA )
!
! 時間積分を行い,
! 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から
! 時刻 $ t+\Delta t $ の物理量を計算します.
!
! 時間積分法にはリープフロッグスキームを用いています.
! ただし拡散項には前方差分を用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! With time integration, calculate physical values at $ t+\Delta t $
! from tendency at $ t $ and physical values at $ t-\Delta t $ .
!
! Leap-frog scheme is used as time integration scheme.
! And forward difference is used to diffusion terms.
! By default, semi-implicit scheme is applied to gravitational terms
! for extension of $ \Delta t $ .
! Explicit scheme can be applied to gravitational terms by changing
! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 座標データ設定
! Axes data settings
!
use axesset, only: z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! LU 分解法により連立 1 次方程式を解くための関数
! Functions to solve linear simultaneous equation by LU decomposition
!
use lumatrix, only: LUSolve
! 文字列操作
! Character handling
!
use dc_string, only: LChar
implicit none
real(DP), intent(in):: wz_DVorDtN ((nmax+1)**2, 1:kmax)
! $ \DD{\zeta}{t} (t) $ . 渦度変化 (スペクトル).
! Vorticity tendency (spectral)
real(DP), intent(in):: wz_DDivDtN ((nmax+1)**2, 1:kmax)
! $ \DD{D}{t} (t) $ . 発散変化 (スペクトル).
! Divergence tendency (spectral)
real(DP), intent(in):: wz_DTempDtN ((nmax+1)**2, 1:kmax)
! $ \DD{T}{t} (t) $ . 温度変化 (スペクトル).
! Temperature tendency (spectral)
real(DP), intent(in):: wz_DQVapDtN ((nmax+1)**2, 1:kmax)
! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP), intent(in):: w_DPiDtN ((nmax+1)**2)
! $ \DD{p_s}{t} (t) $ . 地表面気圧変化 (スペクトル).
! Surface pressure tendency (spectral)
real(DP), intent(in):: wz_VorB ((nmax+1)**2, 1:kmax)
! $ \zeta (t-\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(in):: wz_DivB ((nmax+1)**2, 1:kmax)
! $ D (t-\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(in):: wz_TempB ((nmax+1)**2, 1:kmax)
! $ T (t-\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(in):: w_PiB ((nmax+1)**2)
! $ \pi = \ln p_s (t-\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP), intent(in):: wz_QVapB ((nmax+1)**2, 1:kmax)
! $ q (t-\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP), intent(out):: wz_VorA ((nmax+1)**2, 1:kmax)
! $ \zeta (t+\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(out):: wz_DivA ((nmax+1)**2, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(out):: wz_TempA ((nmax+1)**2, 1:kmax)
! $ T (t+\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(out):: w_PiA ((nmax+1)**2)
! $ \pi = \ln p_s (t+\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP), intent(out):: wz_QVapA ((nmax+1)**2, 1:kmax)
! $ q (t+\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
! 作業変数
! Work variables
!
real(DP):: wz_siTempW ((nmax+1)**2, 1:kmax)
! 温度 (セミインプリシット法のための作業変数).
! Temperature (work variable for semi-implicit scheme)
real(DP):: wz_siDTempDtW ((nmax+1)**2, 1:kmax)
! $ \DD{T}{t} (t) $ . 温度変化 (スペクトル) の作業変数.
! Temperature tendency (spectral) work variable
real(DP):: w_siPiW ((nmax+1)**2)
! $ \pi $ (セミインプリシット法のための作業変数).
! $ \pi $ (work variable for semi-implicit scheme)
real(DP):: w_siDPiDtW ((nmax+1)**2)
! $ \DD{p_s}{t} (t) $ . 地表面気圧変化 (スペクトル).
! Surface pressure tendency (spectral)
real(DP):: wz_siPhiW ((nmax+1)**2, 1:kmax)
! $ \Phi = \underline{W} \overline{ \Dvect{T} }^{t}$ .
! (セミインプリシット法のための作業変数).
! (Work variable for semi-implicit scheme)
real(DP):: wz_siDivAvrTime ((nmax+1)**2, 1:kmax)
! 時間平均の $ \Dvect{D} $ (セミインプリシット法のための作業変数).
! Time average $ \Dvect{D} $ (work variable for semi-implicit scheme)
integer:: k, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
! 非重力波項の仮積分
! Integration non gravitational terms temporarily
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
wz_siTempW = wz_TempB * ( 1.0_DP - DelTime * wz_DiffTherm ) + wz_DTempDtN * DelTime
w_siPiW = w_PiB + w_DPiDtN * DelTime
end select
! ジオポテンシャルの計算
! Calculate geopotential
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
wz_siPhiW (:,1) = CpDry * z_HydroAlpha(1) * wz_siTempW (:,1)
do k = 2, kmax
wz_siPhiW (:,k) = wz_siPhiW(:,k-1) + CpDry * z_HydroAlpha(k) * wz_siTempW (:,k) + CpDry * z_HydroBeta (k-1) * wz_siTempW (:,k-1)
end do
case ('explicit')
end select
! 発散方程式の右辺の計算
! Calculate right side of divergence equation
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
do k = 1, kmax
wz_siDivAvrTime(:,k) = wz_DivB(:,k) * ( 1.0_DP + DelTime * wz_DiffVorDiv(:,k) ) * ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(:,k) ) + DelTime * ( wz_DDivDtN(:,k) - wz_rn(:,k) / RPlanet**2 * ( wz_siPhiW(:,k) + ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(:,k) ) * ( w_Phis + w_siPiW * z_siMtxG(k) ) ) )
end do
end select
! 時間平均の $ \Dvect{D} $ を LU 行列で解く
! Solve time average $ \Dvect{D} $ with LU matrix
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
wz_siDivAvrTime = LUSolve( wzz_siMtxLU, wz_siMtxPiv, wz_siDivAvrTime )
end select
! 温度, 地表気圧の時間変化項の計算
! Calculate tendency of temperature and surface pressure
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
wz_siDTempDtW = wz_DTempDtN
do k = 1, kmax
do kk = 1, kmax
wz_siDTempDtW(:,k) = wz_siDTempDtW(:,k) - zz_siMtxH(k,kk) * wz_siDivAvrTime(:,kk)
end do
end do
w_siDPiDtW = w_DPiDtN
do kk = 1, kmax
w_siDPiDtW = w_siDPiDtW - z_DelSigma(kk) * wz_siDivAvrTime(:,kk)
end do
end select
! 時間積分. 拡散は前方差分
! Time integration. Forward difference is applied to diffusion
!
! 渦度 ; Vorticity
!
wz_VorA = ( wz_VorB + wz_DVorDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffVorDiv )
! 発散 ; Divergence
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
wz_DivA = wz_siDivAvrTime * 2.0_DP - wz_DivB
case ('explicit')
wz_DivA = ( wz_DivB + wz_DDivDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffVorDiv )
end select
! 温度 ; Temperature
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
wz_TempA = ( wz_TempB + wz_siDTempDtW * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm )
case ('explicit')
wz_TempA = ( wz_TempB + wz_DTempDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm )
end select
! 比湿 ; Specific humidity
!
wz_QVapA = ( wz_QVapB + wz_DQVapDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm )
! 地表面気圧 ; Surface pressure
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
w_PiA = w_PiB + w_siDPiDtW * DelTime * 2.0_DP
case ('explicit')
w_PiA = w_PiB + w_DPiDtN * DelTime * 2.0_DP
end select
end subroutine TimeIntegration
| Constant : | |||
| module_name = ‘dynamics_hspl_vas83‘ : | character(*), parameter
|
| Variable : | |||
| r_TempAvrXY(:) : | real(DP), allocatable
|
| Constant : | |||
| version = ’$Name: dcpam5-20081129 $’ // ’$Id: dynamics_hspl_vas83.F90,v 1.1 2008-10-09 10:15:41 morikawa Exp $’ : | character(*), parameter
|
| Variable : | |||
| wz_DiffTherm(:,:) : | real(DP), allocatable
|
| Variable : | |||
| wz_DiffVorDiv(:,:) : | real(DP), allocatable
|
| Variable : | |||
| wz_rn(:,:) : | real(DP), allocatable
|
| Variable : | |||
| wz_siMtxPiv(:,:) : | integer, allocatable
|
| Variable : | |||
| wzz_siMtxLU(:,:,:) : | real(DP), allocatable
|
| Variable : | |||
| xy_Cori(:,:) : | real(DP), allocatable
|
| Variable : | |||
| z_HydroAlpha(:) : | real(DP), allocatable
|
| Variable : | |||
| z_HydroBeta(:) : | real(DP), allocatable
|
| Variable : | |||
| z_TempAvrXY(:) : | real(DP), allocatable
|
| Variable : | |||
| z_TempInpolA(:) : | real(DP), allocatable
|
| Variable : | |||
| z_TempInpolB(:) : | real(DP), allocatable
|
| Variable : | |||
| z_TempInpolKappa(:) : | real(DP), allocatable
|
| Variable : | |||
| zz_siMtxGCt(:,:) : | real(DP), allocatable
|