!= 力学過程 (スペクトル法, Arakawa and Suarez (1983))
!
!= Dynamical process (Spectral method, Arakawa and Suarez (1983))
!
! Authors:: Yasuhiro MORIKAWA, Satoshi NODA, Yoshiyuki O. TAKAHASHI, Shin-ichi Takehiro
! Version:: $Id: dynamics_hspl_vas83.F90,v 1.83 2026/03/19 20:00:00 takepiro Exp $
! Tag Name:: $Name: $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License:: See COPYRIGHT[link:../../../COPYRIGHT]
!
module dynamics_hspl_vas83
!
!= 力学過程 (スペクトル法, Arakawa and Suarez (1983))
!
!= Dynamical Core (Spectral method, Arakawa and Suarez (1983))
!
! Note that Japanese and English are described in parallel.
!
! 力学過程を演算するモジュールです.
! 水平離散化にスペクトル法 (Bourke, 1988) を,
! 鉛直離散化には Arakawa and Suarez (1983) を用いています.
!
! 時間積分法にはリープフロッグスキームを用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! This is a dynamical core module. Spectral method (Bourke, 1988)
! (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
! in order to enlarge the value of $ \Delta t $ .
! Explicit scheme can be applied to gravitational terms by changing
! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
!
!== Procedures List
!
! DynamicsHSplVAS83 :: 力学計算
! DynamicsHSplVAS83Init :: 初期化
! DynamicsHSplVAS83Finalize :: 終了処理 (モジュール内部の変数の割り付け解除)
! --------------------- :: ------------
! DynamicsHSplVAS83 :: Calculate dynamics
! DynamicsHSplVAS83Init :: Initialization
! DynamicsHSplVAS83Finalize :: Termination (deallocate variables in this module)
!
!== NAMELIST
!
! NAMELIST#dynamics_hspl_vas83_nml
!
!== References
!
! * Bourke, W. P., 1988:
! Spectral methods in global climate and weather
! prediction models.
! Physically Based Modelling and Simulation of
! Climates and Climatic Change. Part I.,
! M.E. Schlesinger (ed.),
! Kluwer Academic Publishers, Dordrecht, 169--220.
!
! * Arakawa, A., Suarez, M. J., 1983:
! Vertical differencing of the primitive equations
! in sigma coordinates.
! Mon. Wea. Rev., 111, 34--35.
! モジュール引用 ; USE statements
!
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: DP, & ! 倍精度実数型. Double precision.
& TOKEN ! キーワード. Keywords.
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 宣言文 ; Declaration statements
!
implicit none
private
! 公開手続き
! Public procedure
!
public :: DynamicsHSplVAS83
public :: DynamicsHSplVAS83Init
public :: DynamicsHSplVAS83Finalize
! 公開変数
! Public variables
!
logical, save, public:: wa_module_inited = .false.
! 非公開変数
! Private variables
!
logical, save :: dynamics_hspl_vas83_inited = .false.
! 初期設定フラグ.
! Initialization flag
integer , save :: IDTimeIntegScheme
! ID for used time integration scheme
integer , parameter :: IDTimeIntegSchemeExplicit = 0
integer , parameter :: IDTimeIntegSchemeSemiIMplicit = 1
real(DP), save:: DelTimeSave
! 前回の $ \Delta t $ [s].
! 陰解法のための係数の再生成の必要性の
! チェックに使用する.
!
! $ \Delta t $ [s] at previous step
! This is used to check necessity of
! regeneration of coefficients for
! implicit method.
! 水平拡散係数
! Coefficient of horizontal diffusion
!
real(DP), save, allocatable:: w_LaplaEigVal (:)
! $ \nabla^2_{\sigma} = -n (n+1) $ .
! ラプラシアンの固有値.
! Eigen value of laplacian operator.
real(DP), save, allocatable:: wz_DisCoefM (:,:)
! 運動量の水平拡散係数とスポンジ層減衰係数の和
! Damping coefficient for momentum by horizonatl diffusion
! and a sponge layer
real(DP), save, allocatable:: wz_DisCoefH (:,:)
! 熱の水平拡散係数とスポンジ層減衰係数の和
! Damping coefficient for heat by horizonatl diffusion
! and a sponge layer
real(DP), save, allocatable:: w_DisCoefQ (:)
! 成分の水平拡散係数
! Damping coefficient for composition by horizonatl diffusion
! NonLinearOnGrid 等で使用する係数
! Coefficients for "NonLinearOnGrid", etc.
!
real(DP), allocatable:: xy_SinLat (:,:)
! $ \sin \varphi $ .
real(DP), allocatable:: xy_CosLat (:,:)
! $ \cos \varphi $ .
real(DP), allocatable:: xy_Cori (:,:)
! $ f\equiv 2\Omega\sin\varphi $ .
! コリオリパラメータ. Coriolis parameter
real(DP), allocatable:: z_HydroAlpha (:)
! $ \alpha $ . 静水圧の式の係数.
! Coefficient in hydrostatic equation.
real(DP), allocatable:: z_HydroBeta (:)
! $ \beta $ . 静水圧の式の係数.
! Coefficient in hydrostatic equation.
real(DP), allocatable:: z_TInpCoefA (:)
! $ a $ . 温度鉛直補間の係数.
! Coefficient for vertical interpolation of temperature
real(DP), allocatable:: z_TInpCoefB (:)
! $ b $ . 温度鉛直補間の係数.
! Coefficient for vertical interpolation of temperature
real(DP), allocatable:: z_TInpCoefK (:)
! $ \kappa $ . 温度鉛直補間の係数.
! Coefficient for vertical interpolation of temperature
real(DP), allocatable:: z_RefTemp (:)
! $ \overline{T} $ . 基準温度 (整数レベル).
! Reference temperature on full sigma level
real(DP), allocatable:: r_RefTemp (:)
! $ \hat{\overline{T}} $ . 基準温度 (半整数レベル).
! Reference temperature on half sigma level
! TimeIntegration で使用する係数
! Coefficients for "TimeIntegration"
!
real(DP), save, allocatable:: z_siMtxC (:)
! $ C = \sigma $
real(DP), save, allocatable:: z_siMtxG (:)
! $ G = C_p \kappa T $
real(DP), save, allocatable:: zz_siMtxH (:,:)
! $ h = QS - R $
real(DP), save, allocatable:: wz_siMtxDi (:,:)
! $ D = ( 1 - 2 \Delta t D_h )^{-1}$
! NOTE: This is a diagonal matrix.
real(DP), save, allocatable:: zz_siMtxDiH(:,:)
! $ Di h $
real(DP), save, allocatable:: wzz_siMtxWDiH(:,:,:)
! $ W Di h $
real(DP), save, allocatable:: zz_siMtxGCt (:,:)
! $ G C^{T} $ ( $ C = \Delta \sigma $ )
real(DP), save, allocatable:: wzz_siMtxLU(:,:,:)
! セミインプリシット行列の LU 分解.
! LU decomposition for semi-implicit matrix
integer , save, allocatable:: wz_siMtxPiv(:,:)
! セミインプリシット行列のピボット.
! Pivot for semi-implicit matrix
! サブルーチン SemiImplMatrix 中で使用する係数
! Coefficients used in a subroutine "SemiImplMatrix"
!
real(DP), save, allocatable:: zz_siMtxW (:,:)
! $ W $ .
! 発散の式での線形重力波項の効果による温度補正の係数.
! Coefficient for correction of temperature
! by effort of linear gravitational terms
real(DP), save, allocatable:: zz_siMtxQ (:,:)
! $ Q = \DD{T}{\sigma} $
real(DP), save, allocatable:: zz_siMtxS (:,:)
! $ S = \DD{\sigma}{t} $
real(DP), save, 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, save, allocatable:: nmo (:,:,:)
!!$ ! スペクトルの添字順番.
!!$ ! Spectral subscript expression
!!$ real(DP), save, allocatable:: wzz_siMtxM (:,:,:)
!!$ ! 行列 $ \underline{M} $.
!!$ ! Matrix $ \underline{M} $
!!$ integer, save, allocatable:: z_siMtxPivWork(:)
!!$ ! 行列のピボット作成の作業変数.
!!$ ! Work variable for pivot
logical , save :: FlagDivDamp
! Flag for divergence damping
real(DP), save :: DivDampPeriod
! Period for divergence damping application in unit of
! second
real(DP), save :: DivDampTime0
! Time for starting divergence damping in unit of second
logical , save :: FlagMassHorDifCor
! Flag for correction of horizontal diffusion of mass
logical , save :: FlagTotMassFix
! Flag for fixing total mass
logical , save :: FlagSLTT
!!$ integer , save :: IDTTMethod
!!$ integer , parameter :: IDTTMethodSp = 1
!!$ integer , parameter :: IDTTMethodSL = 2
!!$ ! セミラグランジュ法による物質移流フラグ。
!!$ ! Flag for using Semi-Lagrangian method for tracer transport.
!!$ integer , parameter :: IDTTMethodFV = 3
logical , save :: FlagAdvTest
!
! Flag for advection test.
! When this is true, U, V, T, and Ps are not
! calculated prognostically.
character(*), parameter:: module_name = 'dynamics_hspl_vas83'
! モジュールの名称.
! Module name
character(*), parameter:: version = &
& '$Name: $' // &
& '$Id: dynamics_hspl_vas83.F90,v 1.84 2026/03/19 20:00:00 takepiro Exp $'
! モジュールのバージョン
! Module version
contains
!--------------------------------------------------------------------------------------
subroutine DynamicsHSplVAS83( &
& xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, & ! (in)
& xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, & ! (in)
& xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, & ! (in)
& xy_SurfHeight, & ! (in)
& xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA, & ! (out)
& xyz_ArgOMG & ! (out)
& )
!
! 力学過程の演算を行い, 与えられた $ t-\Delta t $ および $ t $ の
! 東西風速 (xyz_UB, xyz_UN), 南北風速 (xyz_VB, xyz_VN),
! 温度 (xyz_TempB, xyz_TempN), 質量混合比 (xyzf_QMixB, xyzf_QMixN),
! 地表面気圧 (xyz_PsB, xyz_PsN), および地表面高度 (xy_SurfHeight) から,
! $ t+\Delta t $ の
! 東西風速 (xyz_UA), 南北風速 (xyz_VA), 温度 (xyz_TempA),
! 質量混合比 (xyzf_QMixA), 地表面気圧 (xyz_PsA) を返します.
!
! 別の物理プロセスによる渦度, 発散, 温度, 比湿の変化を,
! 力学過程の変化に足して次のステップを計算する場合には,
! それらの変化を xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy
! に与えてください.
!
! 時間積分法にはリープフロッグスキームを用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! Calculate dynamical processes.
! Input data are
! Eastward wind (xyz_UB, xyz_UN),
! northward wind (xyz_VB, xyz_VN),
! temperature (xyz_TempB, xyz_TempN),
! mass mixing ratio (xyzf_QMixB, xyzf_QMixN),
! surface pressure (xyz_PsB, xyz_PsN) at $ t-\Delta t $ and $ t $,
! and surface height (xy_SurfHeight).
! Eastward wind (xyz_UA), northward wind (xyz_VA),
! temperature (xyz_TempA),
! mass mixing ratio (xyzf_QMixA), surface pressure (xyz_PsA)
! are returned.
!
! In order to add tendencies of vorticity, divergence,
! temperature and specific humidity calculated by
! physical processes other than those by
! this dynamical process,
! these tendencies should be given to
! "xyz_DUDtPhy", "xyz_DVDtPhy", "xyz_DTempDtPhy", "xyzf_DQMixDtPhy"
!
! 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, &
! $ a $ [m].
! 惑星半径.
! Radius of planet
& CpDry, &
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
& Grav
! $ g $ [m s-2].
! 重力加速度.
! Gravitational acceleration
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax, &
! 成分の数
! Number of composition
& IndexH2OVap, &
& CompositionInqFlagAdv
! 時刻管理
! Time control
!
use timeset, only: &
& DelTime, & ! $ \Delta t $ [s]
& TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $.
& TimesetClockStart, TimesetClockStop
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& r_Sigma, &
! $ \sigma $ レベル (半整数).
! Half $ \sigma $ level
& z_Sigma ! $ \sigma $ レベル (整数).
! Full $ \sigma $ level
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef AXISYMMETRY
use wa_zonal_module, only: &
& wa_Lapla_wa, &
& w_Lapla_w, &
& wa_LaplaInv_wa, &
& wa_DivLambda_xya, &
& wa_DivMu_xya, &
& xy_GradLambda_w, &
& xya_GradLambda_wa, &
& xy_GradMu_w, &
& xya_GradMu_wa, &
& w_xy, xy_w, &
& wa_xya, xya_wa
#elif LIB_MPI
! modify for bug of Intel fortran (2025/09/11 S.Takehiro)
! use ua_mpi_module, only: &
! & wa_Lapla_wa => ua_Lapla_ua, &
! & w_Lapla_w => u_Lapla_u, &
! & wa_LaplaInv_wa => ua_LaplaInv_ua, &
! & wa_DivLambda_xya => ua_DivLambda_pva, &
! & wa_DivMu_xya => ua_DivMu_pva, &
! & xy_GradLambda_w => pv_GradLambda_u, &
! & xya_GradLambda_wa => pva_GradLambda_ua, &
! & xy_GradMu_w => pv_GradMu_u, &
! & xya_GradMu_wa => pva_GradMu_ua, &
! & w_xy => u_pv, &
! & xy_w => pv_u, &
! & wa_xya => ua_pva, &
! & xya_wa => pva_ua
use ua_mpi_module_base, only: &
& w_xy => u_pv, &
& xy_w => pv_u, &
& wa_xya => ua_pva, &
& xya_wa => pva_ua
use ua_mpi_module_deriv, only: &
& wa_Lapla_wa => ua_Lapla_ua, &
& w_Lapla_w => u_Lapla_u, &
& wa_LaplaInv_wa => ua_LaplaInv_ua, &
& wa_DivLambda_xya => ua_DivLambda_pva, &
& wa_DivMu_xya => ua_DivMu_pva, &
& xy_GradLambda_w => pv_GradLambda_u, &
& xya_GradLambda_wa => pva_GradLambda_ua, &
& xy_GradMu_w => pv_GradMu_u, &
& xya_GradMu_wa => pva_GradMu_ua
#else
use wa_module, only: &
& wa_Lapla_wa, &
& w_Lapla_w, &
& wa_LaplaInv_wa, &
& wa_DivLambda_xya, &
& wa_DivMu_xya, &
& xy_GradLambda_w, &
& xya_GradLambda_wa, &
& xy_GradMu_w, &
& xya_GradMu_wa, &
& 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.
! 積分と平均の操作
! Operation for integral and average
!
use intavr_operate, only: IntLonLat_xy
! 質量の補正
! Mass fixer
!
use mass_fixer, only: MassFixer, MassFixerColumn
! セミラグランジュ法による物質移流計算
! Semi-Lagrangian method for tracer transport
!
use sltt, only: SLTTMain
!!$ ! 物質移流 (有限体積法)
!!$ ! Tracer Transport (finite volume method)
!!$ !
!!$ use tt_fv, only: TTFVMain
!
! Utility module for advection test
!
use adv_test, only: AdvTestSetHorVels, AdvTestSetVerVel
!
! 水平スペクトル解析
! horizontal spectral analysis
!
use dynamics_diag_spectrum, only : DynamicsDiagSpAnalysis
! 宣言文 ; 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):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ 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):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ 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_DUDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{u}{t}\right)^{phy} $ .
! 外力項 (物理過程) による東西風速変化.
! Eastward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DVDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{v}{t}\right)^{phy} $ .
! 外力項 (物理過程) による南北風速変化.
! Northward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DTempDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{T}{t}\right)^{phy} $ .
! 外力項 (物理過程) による温度変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \left(\DP{q}{t}\right)^{phy} $ .
! 外力項 (物理過程) による比湿変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xy_SurfHeight (0:imax-1, 1:jmax)
! $ z_s $ . 地表面高度.
! Surface height.
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):: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿. Specific humidity
real(DP), intent(out):: xy_PsA (0:imax-1, 1:jmax)
! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure
real(DP), intent(out), optional :: xyz_ArgOMG(0:imax-1, 1:jmax, 1:kmax)
!
! OMEGA, DP/Dt
! 作業変数
! Work variables
!
real(DP):: xyz_UCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ U (t) = u (t) \cos \varphi $ .
real(DP):: xyz_VCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ V (t) = v (t) \cos \varphi $ .
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):: wz_DivN (lmax, 1:kmax)
! $ D (t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempN (lmax, 1:kmax)
! $ T (t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiN (lmax)
! $ \pi $ スペクトル
real(DP):: xy_GradLambdaPiN (0:imax-1, 1:jmax)
! $ \DP{\pi}{\lambda} $
real(DP):: xy_GradMuPiN (0:imax-1, 1:jmax)
! $ (1-\mu^2) \DP{\pi}{\mu} $
real(DP):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{v} \cdot \nabla \pi $ .
! $ \pi $ の移流. Advection of $ \pi $
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_TempNonLinearN (0:imax-1, 1:jmax, 1:kmax)
! $ H (t) $ . 温度時間変化項.
! Temperature tendency
real(DP):: xyzf_QMixNonLinearN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ 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_SigDotN (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流. Vertical flow
real(DP):: xy_DPiDtNG (0:imax-1, 1:jmax)
! $ Z $ . 地表面気圧時間変化非重力波項.
! Non-gravity wave component of surface pressure tendency
real(DP):: xyzf_QMixUAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Uq (t) $ . 比湿東西移流項.
! Eastward advection of specific humidity
real(DP):: xyzf_QMixVAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Vq (t) $ . 比湿南北移流項.
! Northward advection of specific humidity
real(DP):: wz_DVorDtNG (lmax, 1:kmax)
! $ \left( \DD{\zeta}{t} (t) \right)^{NG}$ . 渦度変化の非重力波成分 (スペクトル).
! Non-gravity wave component of vorticity tendency (spectral)
real(DP):: wz_DDivDtNG (lmax, 1:kmax)
! $ \left( \DD{D}{t} (t) \right)^{NG}$ . 発散変化の非重力波成分 (スペクトル).
! Non-gravity wave component of divergence tendency (spectral)
real(DP):: wz_DTempDtNG (lmax, 1:kmax)
! $ \left( \DD{T}{t} (t) \right)^{NG}$ . 温度変化の非重力波成分 (スペクトル).
! Non-gravity wave component of temperature tendency (spectral)
real(DP):: wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax)
! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP):: w_DPiDtNG(lmax)
! $ \left( \DD{p_s}{t} (t) \right)^{NG}$ . 地表面気圧変化費重力波項 (スペクトル).
! Non-gravity wave component of surface pressure tendency (spectral)
real(DP):: xyz_UCosLatB (0:imax-1, 1:jmax, 1:kmax)
! $ U (t-\Delta t) = u (t-\Delta t) \cos \varphi $ .
real(DP):: xyz_VCosLatB (0:imax-1, 1:jmax, 1:kmax)
! $ V (t-\Delta t) = v (t-\Delta t) \cos \varphi $ .
real(DP):: wz_VorB (lmax, 1:kmax)
! $ \zeta (t-\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivB (lmax, 1:kmax)
! $ D (t-\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempB (lmax, 1:kmax)
! $ T (t-\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiB (lmax)
! $ \pi = \ln p_s (t-\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wzf_QMixB(lmax, 1:kmax, 1:ncmax)
! $ q (t-\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP):: xyz_DUDtPhyCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{u}{t}\right)^{phy} \cos \varphi $ .
real(DP):: xyz_DVDtPhyCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{v}{t}\right)^{phy} \cos \varphi $ .
real(DP):: wz_VorA (lmax, 1:kmax)
! $ \zeta (t+\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivA (lmax, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempA (lmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiA (lmax)
! $ \pi = \ln p_s (t+\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wzf_QMixA(lmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP):: wz_Psi (lmax, 1:kmax)
! $ \psi $ . 流線関数. Streamline function
real(DP):: wz_Chi (lmax, 1:kmax)
! $ \chi $ . ポテンシャル. Potential
real(DP):: xyz_UCosLatA (0:imax-1, 1:jmax, 1:kmax)
! $ U (t+\Delta t) = u (t+\Delta t) \cos \varphi $ .
real(DP):: xyz_VCosLatA (0:imax-1, 1:jmax, 1:kmax)
! $ V (t+\Delta t) = v (t+\Delta t) \cos \varphi $ .
real(DP):: wz_VorDiffA (lmax, 1:kmax)
! $ \mathscr{D}(\zeta) $ .
! 運動量水平拡散による渦度変化 (スペクトル).
! Vorticity tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_DivDiffA (lmax, 1:kmax)
! $ \mathscr{D}(D) $ .
! 運動量水平拡散による発散変化 (スペクトル).
! Divergence tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_PsiDiff (lmax, 1:kmax)
! 運動量水平拡散による流線関数 $ \psi $ 変化
! Streamline function tendency by
! horizontal momentum diffusion
real(DP):: wz_ChiDiff (lmax, 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
! 地表面高度関連
! Surface height etc.
!
real(DP):: w_SurfGeoPot (lmax)
! $ \Phi_s $ . 地表ジオポテンシャル.
! Surface geo-potential
! Variables for mass fixing
!
real(DP):: MassB
real(DP):: MassA
real(DP):: xyr_PressA(0:imax-1, 1:jmax, 0:kmax)
real(DP):: xyr_PressB(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyz_PressA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_HDifPsA (0:imax-1, 1:jmax)
real(DP) :: xyz_HDifPressA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyzf_DQMixDtHDCor(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
integer :: kp, kn
real(DP) :: xyz_OMG (0:imax-1, 1:jmax, 1:kmax)
!
! OMEGA, DP/Dt
real(DP) :: xyzf_QMixTentative(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
logical :: FlagExistNonAdvComp
real(DP) :: xyz_UTracerAdv(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_VTracerAdv(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyr_AdvTestStreamFunc(0:imax-1, 1:jmax, 0:kmax)
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. dynamics_hspl_vas83_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
call SemiImplMatrix
! 格子点値をスペクトル値へ ( $ t-\Delta t$ )
! Exchange grid values to spectral values ( $ t-\Delta t$ )
!
! 渦度発散の計算
! Calculate vorticity and divergence
do k = 1, kmax
xyz_UCosLatB(:,:,k) = xyz_UB(:,:,k) * xy_CosLat
xyz_VCosLatB(:,:,k) = xyz_VB(:,:,k) * xy_CosLat
end do
wz_VorB = ( wa_DivLambda_xya( xyz_VCosLatB ) &
& - wa_DivMu_xya ( xyz_UCosLatB ) ) / RPlanet
wz_DivB = ( wa_DivLambda_xya( xyz_UCosLatB ) &
& + wa_DivMu_xya ( xyz_VCosLatB ) ) / RPlanet
!
wz_TempB = wa_xya( xyz_TempB )
!
w_PiB = w_xy( log( xy_PsB ) )
! 格子点値をスペクトル値へ ( $ t $ )
! Exchange grid values to spectral values ( $ t $ )
!
! 渦度発散の計算
! Calculate vorticity and divergence
!
do k = 1, kmax
xyz_UCosLatN(:,:,k) = xyz_UN(:,:,k) * xy_CosLat
xyz_VCosLatN(:,:,k) = xyz_VN(:,:,k) * xy_CosLat
end do
xyz_VorN = xya_wa( wa_DivLambda_xya( xyz_VCosLatN ) &
& - wa_DivMu_xya ( xyz_UCosLatN ) ) / RPlanet
wz_DivN = ( wa_DivLambda_xya( xyz_UCosLatN ) &
& + wa_DivMu_xya ( xyz_VCosLatN ) ) / RPlanet
xyz_DivN = xya_wa( wz_DivN )
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeExplicit )
wz_TempN = wa_xya( xyz_TempN )
end select
!
xy_PiN = log( xy_PsN )
w_PiN = w_xy( xy_PiN )
xy_GradLambdaPiN = xy_GradLambda_w( w_PiN ) / RPlanet
xy_GradMuPiN = xy_GradMu_w ( w_PiN ) / RPlanet
! 地表ジオポテンシャルの計算
! Calculate surface geo-potential
!
w_SurfGeoPot = w_xy( Grav * xy_SurfHeight )
! 格子点上での非線形力学項の計算
! Calculate non-linear dynamical terms on grid points
!
call NonLinearOnGrid( &
& xyz_UCosLatN, xyz_VCosLatN, & ! (in)
& xyz_VorN, xyz_DivN, & ! (in)
& xyz_TempN, xyzf_QMixN(:,:,:,IndexH2OVap), & ! (in)
& xy_GradLambdaPiN, xy_GradMuPiN, & ! (in)
!
& xyz_PiAdv, & ! (out)
& xyz_UAdvN, xyz_VAdvN, & ! (out)
& xyz_TempNonLinearN, & ! (out)
& xyz_KinEngyN, & ! (out)
& xyz_TempUAdvN, xyz_TempVAdvN, & ! (out)
& xyr_SigDotN, xy_DPiDtNG & ! (out)
& )
! 非線形項と物理過程による時間変化項の和をスペクトル変換
! Spectral transformation of sum of non-linear term and tendencies by physical
! processes
!
do k = 1, kmax
xyz_DUDtPhyCosLat(:,:,k) = xyz_DUDtPhy(:,:,k) * xy_CosLat
xyz_DVDtPhyCosLat(:,:,k) = xyz_DVDtPhy(:,:,k) * xy_CosLat
end do
wz_DVorDtNG = ( wa_DivLambda_xya( xyz_VAdvN + xyz_DVDtPhyCosLat ) &
& - wa_DivMu_xya ( xyz_UAdvN + xyz_DUDtPhyCosLat ) ) &
& / RPlanet
wz_DDivDtNG = ( wa_DivLambda_xya( xyz_UAdvN + xyz_DUDtPhyCosLat ) &
& + wa_DivMu_xya ( xyz_VAdvN + xyz_DVDtPhyCosLat ) ) &
& / RPlanet &
& - wa_Lapla_wa( wa_xya(xyz_KinEngyN) ) / RPlanet**2
wz_DTempDtNG = - ( wa_DivLambda_xya( xyz_TempUAdvN ) &
& + wa_DivMu_xya ( xyz_TempVAdvN ) ) &
& / RPlanet &
& + wa_xya( xyz_TempNonLinearN + xyz_DTempDtPhy )
w_DPiDtNG = w_xy( xy_DPiDtNG )
! 時間積分
! Time integration
!
call TimeIntegration( &
& w_SurfGeoPot, & ! (in)
& wz_DVorDtNG, wz_DDivDtNG, wz_DTempDtNG, w_DPiDtNG, & ! (in)
& wz_DivN, wz_TempN, w_PiN, & ! (in)
& wz_VorB, wz_DivB, wz_TempB, w_PiB, & ! (in)
& wz_VorA, wz_DivA, wz_TempA, w_PiA & ! (out)
& )
! Divergence damping
!
call DivergenceDamping( &
& wz_DivA & ! (inout)
& )
! スペクトル値を格子点値へ ( $ 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_UCosLatA = &
& ( xya_GradLambda_wa( wz_Chi ) &
& - xya_GradMu_wa ( wz_Psi ) ) / RPlanet
xyz_VCosLatA = &
& ( xya_GradLambda_wa( wz_Psi ) &
& + xya_GradMu_wa ( wz_Chi ) ) / RPlanet
do k = 1, kmax
xyz_UA(:,:,k) = xyz_UCosLatA(:,:,k) / xy_CosLat
xyz_VA(:,:,k) = xyz_VCosLatA(:,:,k) / xy_CosLat
end do
xyz_TempA = xya_wa( wz_TempA )
xy_PsA = exp( xy_w( w_PiA ) )
! 水平拡散とスポンジ層による散逸の補正
! Correction for dissipation of kinetic enery by horizontal diffusion and sponge
! layer
!
! 水平拡散とスポンジ層による渦度発散の時間変化
! Vorticity and divergence tendency by horizontal diffusion and damping in a sponge
! layer
!
wz_VorDiffA = wz_VorA * wz_DisCoefM
wz_DivDiffA = wz_DivA * wz_DisCoefM
! 水平拡散とスポンジ層による散逸の摩擦熱補正
! Correction for internal energy by horizontal diffusion and damping in a sponge
! layer
!
wz_PsiDiff = wa_LaplaInv_wa( wz_VorDiffA ) * RPlanet**2
wz_ChiDiff = wa_LaplaInv_wa( wz_DivDiffA ) * RPlanet**2
xyz_UDiff = &
& ( xya_GradLambda_wa( wz_ChiDiff ) &
& - xya_GradMu_wa ( wz_PsiDiff ) ) / RPlanet
xyz_VDiff = &
& ( xya_GradLambda_wa( wz_PsiDiff ) &
& + xya_GradMu_wa ( wz_ChiDiff ) ) / RPlanet
do k = 1, kmax
xyz_TempA(:,:,k) = xyz_TempA(:,:,k) &
& - ( xyz_UA(:,:,k) * xyz_UDiff(:,:,k) &
& + xyz_VA(:,:,k) * xyz_VDiff(:,:,k) ) &
& / xy_CosLat / CpDry * 2.0_DP * DelTime
end do
! Mass fixer
! Total Mass
! NOTICE: It is assumed that the total atmospheric mass does not change.
!
if ( FlagTotMassFix ) then
MassB = IntLonLat_xy( xy_PsB )
MassA = IntLonLat_xy( xy_PsA )
if ( MassA /= 0.0_DP ) then
xy_PsA = MassB / MassA * xy_PsA
end if
end if
! Set velocities for tracer advection
!
xyz_UTracerAdv = xyz_UN
xyz_VTracerAdv = xyz_VN
! 主に移流テストのために使う if 文
!
if ( FlagAdvTest ) then
xyz_UA = xyz_UB
xyz_VA = xyz_VB
xyz_TempA = xyz_TempB
xy_PsA = xy_PsB
xyr_SigDotN = 0.0_DP
!
! Utility module for advection test
!
call AdvTestSetHorVels( &
& xyz_UTracerAdv, xyz_VTracerAdv & ! (out)
& )
! calculation of variables required in spectral advection calculation
do k = 1, kmax
xyz_UCosLatN(:,:,k) = xyz_UTracerAdv(:,:,k) * xy_CosLat
xyz_VCosLatN(:,:,k) = xyz_VTracerAdv(:,:,k) * xy_CosLat
end do
wz_DivN = ( wa_DivLambda_xya( xyz_UCosLatN ) &
& + wa_DivMu_xya ( xyz_VCosLatN ) ) / RPlanet
xyz_DivN = xya_wa( wz_DivN )
call AdvTestSetVerVel( &
& xyr_SigDotN, & ! (out)
& xyr_AdvTestStreamFunc & ! (out) optional
& )
call HistoryAutoPut( TimeN, 'AdvTestStreamFunc', xyr_AdvTestStreamFunc )
end if
if ( FlagSLTT ) then
!!$ select case ( IDTTMethod )
!!$ case ( IDTTMethodSL )
! Semi-Lagrangian advection
do k = 0, kmax
xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
end do
! Mass fixer is included in SLTTMain
call SLTTMain(&
& xyr_PressB, xyr_PressA, &
& xyz_UTracerAdv, xyz_VTracerAdv, xyr_SigDotN, & ! (in)
& xyzf_DQMixDtPhy, & ! (in)
& xyzf_QMixB, & ! (in)
& xyzf_QMixA & ! (out)
&)
!!$ case ( IDTTMethodFV )
!!$ ! Finite volume method
!!$
!!$ do k = 0, kmax
!!$ xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
!!$ xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
!!$ end do
!!$ call TTFVMain( &
!!$ & xyr_PressB, xyr_PressA, &
!!$ & xyz_UTracerAdv, xyz_VTracerAdv, xyr_SigDotN, & ! (in)
!!$ & xyzf_DQMixDtPhy, & ! (in)
!!$ & xyzf_QMixB, & ! (in)
!!$ & xyzf_QMixA & ! (out)
!!$ & )
else
!!$ case ( IDTTMethodSp )
! Spectral advection
do n = 1, ncmax
wzf_QMixB(:,:,n) = wa_xya( xyzf_QMixB(:,:,:,n) )
end do
!
call NonLinearOnGridQMix( &
& xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, & ! (in)
& xyzf_QMixN, & ! (in)
!
& xyzf_QMixNonLinearN, & ! (out)
& xyzf_QMixUAdvN, xyzf_QMixVAdvN & ! (out)
& )
do n = 1, ncmax
wzf_DQMixDtN(:,:,n) = - ( wa_DivLambda_xya( xyzf_QMixUAdvN(:,:,:,n) ) &
& + wa_DivMu_xya ( xyzf_QMixVAdvN(:,:,:,n) ) ) &
& / RPlanet &
& + wa_xya( xyzf_QMixNonLinearN(:,:,:,n) &
& + xyzf_DQMixDtPhy (:,:,:,n) )
end do
call TimeIntegrationQMix( &
& wzf_DQMixDtN, & ! (in)
& wzf_QMixB, & ! (in)
& wzf_QMixA & ! (out)
& )
do n = 1, ncmax
xyzf_QMixA(:,:,:,n) = xya_wa( wzf_QMixA(:,:,n) )
end do
if ( FlagMassHorDifCor ) then
! Correction for horizontal diffusion of constituents.
! This is performed for the aim of correcting the diffusion in the viscinity of
! steep terrain slope.
!
do k = 1, kmax
xyz_PressA(:,:,k) = xy_PsA * z_Sigma(k)
end do
xy_HDifPsA = xy_w( w_DisCoefQ * w_xy( xy_PsA ) )
do k = 1, kmax
xyz_HDifPressA(:,:,k) = z_Sigma(k) * xy_HDifPsA(:,:)
end do
do n = 1, ncmax
do k = 1, kmax
kp = k-1
kn = k+1
if( k == 1 ) kp = 1
if( k == kmax ) kn = kmax
xyzf_DQMixDtHDCor(:,:,k,n) = &
& - ( xyzf_QMixA(:,:,kn,n) - xyzf_QMixA(:,:,kp,n) ) &
& / ( xyz_PressA(:,:,kn) - xyz_PressA(:,:,kp) ) &
& * xyz_HDifPressA(:,:,k)
end do
end do
xyzf_QMixA = xyzf_QMixA + xyzf_DQMixDtHDCor * 2.0_DP * DelTime
endif
! Mass fixer
! Constituents
!
do k = 0, kmax
xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
end do
call MassFixer( &
& xyr_PressA, & ! (in)
& xyzf_QMixA, & ! (inout)
& xyr_PressRef = xyr_PressB, & ! (in) optional
& xyzf_QMixRef = ( xyzf_QMixB+xyzf_DQMixDtPhy*2.0_DP*DelTime ) & ! (in) optional
& )
!!$ end select
end if
! 主に移流テストのために使う if 文
!
!!$ if ( FlagAdvTest ) then
!!$ xyr_SigDotN = 0.0_DP
!!$ end if
! Calculation for non-advection case
! オイラー移流計算のとき、ピーキーな分布の物質を移流させない。計算安定の
! ため。
!
FlagExistNonAdvComp = .false.
do n = 1, ncmax
if ( .not. CompositionInqFlagAdv( n ) ) then
xyzf_QMixA(:,:,:,n) = xyzf_QMixB(:,:,:,n) &
& + xyzf_DQMixDtPhy(:,:,:,n) * 2.0_DP * DelTime
xyzf_QMixTentative(:,:,:,n) = xyzf_QMixA(:,:,:,n)
FlagExistNonAdvComp = .true.
else
xyzf_QMixTentative(:,:,:,n) = 0.0_DP
end if
end do
if ( FlagExistNonAdvComp ) then
! Mass fixer
! Constituents
call MassFixerColumn( &
& xyr_PressB, & ! (in)
& xyzf_QMixTentative, & ! (inout)
& xyr_PressRef = xyr_PressB, & ! (in) optional
& xyzf_QMixRef = ( xyzf_QMixB+xyzf_DQMixDtPhy*2.0_DP*DelTime ) & ! (in) optional
& )
do n = 1, ncmax
if ( .not. CompositionInqFlagAdv( n ) ) then
xyzf_QMixA(:,:,:,n) = xyzf_QMixTentative(:,:,:,n)
end if
end do
end if
! ヒストリデータ出力
! History data output
!
call OutputDiagnosedVariables( &
& xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, & ! (in)
& xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, & ! (in)
& xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA, & ! (in)
& xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, & ! (in)
& xyr_SigDotN, xy_DPiDtNG, & ! (in)
& xyz_PiAdv, & ! (in)
& xyz_OMG & ! (out)
& )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
! 水平スペクトル解析
! Horizontal spectral analysis
!
call DynamicsDiagSpAnalysis( wz_Psi, wz_Chi )
! 計算時間計測再開
! Start measurement of computation time
!
call TimesetClockStart( module_name )
if ( present( xyz_ArgOMG ) ) then
xyz_ArgOMG = xyz_OMG
end if
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine DynamicsHSplVAS83
!--------------------------------------------------------------------------------------
subroutine NonLinearOnGrid( &
& xyz_UCosLatN, xyz_VCosLatN, & ! (in)
& xyz_VorN, xyz_DivN, & ! (in)
& xyz_TempN, xyz_QVapN, & ! (in)
& xy_GradLambdaPiN, xy_GradMuPiN, & ! (in)
!
& xyz_PiAdv, & ! (out)
& xyz_UAdvN, xyz_VAdvN, & ! (out)
& xyz_TempNonLinearN, & ! (out)
& xyz_KinEngyN, & ! (out)
& xyz_TempUAdvN, xyz_TempVAdvN, & ! (out)
& xyr_SigDotN, xy_DPiDtNG & ! (out)
& )
!
! 非線形項 (非重力波項) を格子点上で計算します.
!
! Non-linear terms (non gravitational terms) are calculated on
! grid points
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& r_Sigma, & ! $ \sigma $ レベル (半整数).
! Half $ \sigma $ level
& z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! 物理定数設定
! Physical constants settings
!
use constants, only: &
& CpDry, &
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
& EpsV
! $ \epsilon_v $ .
! 水蒸気分子量比.
! Molecular weight of water vapor
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 格子点設定
! Grid points settings
!
use gridset, only: imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax
! 成分の数
! Number of composition
implicit none
real(DP), intent(in):: xyz_UCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ U (t) = u (t) \cos \varphi $ .
real(DP), intent(in):: xyz_VCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ V (t) = v (t) \cos \varphi $ .
real(DP), intent(in):: xyz_VorN (0:imax-1, 1:jmax, 1:kmax)
! $ \zeta (t) $ . 渦度. Vorticity
real(DP), intent(in):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
! $ D (t) $ . 発散. Divergence
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_GradLambdaPiN (0:imax-1, 1:jmax)
! $ \DP{\pi}{\lambda} $
real(DP), intent(in):: xy_GradMuPiN (0:imax-1, 1:jmax)
! $ (1-\mu^2) \DP{\pi}{\mu} $
real(DP), intent(out):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{v} \cdot \nabla \pi $ .
! $ \pi $ の移流. Advection of $ \pi $
real(DP), intent(out):: xyz_UAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ U_A (t) $ . 東西運動量移流項.
! Eastward advection of momentum
real(DP), intent(out):: xyz_VAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ V_A (t) $ . 南北運動量移流項.
! Northward advection of momentum
real(DP), intent(out):: xyz_TempNonLinearN (0:imax-1, 1:jmax, 1:kmax)
! $ \hat{H} (t) $ . 温度時間変化項.
! Temperature tendency
real(DP), intent(out):: xyz_KinEngyN (0:imax-1, 1:jmax, 1:kmax)
! $ KE (t) $ . 運動エネルギー項.
! Kinetic energy
real(DP), intent(out):: xyz_TempUAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ UT' (t) $ . 温度東西移流項.
! Eastward advection of temperature
real(DP), intent(out):: xyz_TempVAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ VT' (t) $ . 温度南北移流項.
! Northward advection of temperature
real(DP), intent(out):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} (t) $ .
! 鉛直流. Vertical flow
real(DP), intent(out):: xy_DPiDtNG (0:imax-1, 1:jmax)
! $ Z (t) $ . 地表面気圧時間変化非重力波項.
! Non-gravity wave component of surface pressure tendency
!-----------------------------------
! 作業変数
! Work variables
!
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_SigDotNonG (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_VirTemp (0:imax-1, 1:jmax, 1:kmax)
! $ T_v $ .
! 仮温度. Virtual temperature
real(DP):: xyz_VirTempEdd (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 advection of $ \pi $ and divergence
! from the top of the model downward
!
do k = 1, kmax
xyz_PiAdv(:,:,k) = &
& ( xyz_UCosLatN(:,:,k) * xy_GradLambdaPiN &
& + xyz_VCosLatN(:,:,k) * xy_GradMuPiN ) &
& / ( 1.0_DP - xy_SinLat**2 )
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_DivN(:,:,kmax) * z_DelSigma(kmax)
do k = kmax-1, 1, -1
xyz_DivSum(:,:,k) = &
& xyz_DivSum(:,:,k+1) + xyz_DivN(:,:,k) * z_DelSigma(k)
enddo
! 地表面気圧時間変化の非重力波項 $ Z $ の計算
! Calculate non-gravity wave component of surface pressure tendency $ Z $
!
xy_DPiDtNG = - xyz_PiAdvSum(:,:,1)
! $ \dot{\sigma} $ の計算
! Calculate $ \dot{\sigma} $
!
do k = 1, kmax-1
xyr_SigDotN(:,:,k) = &
& r_Sigma(k) * ( xyz_PiAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) &
& - ( xyz_PiAdvSum(:,:,k+1) + xyz_DivSum(:,:,k+1) )
xyr_SigDotNonG(:,:,k) = &
& r_Sigma(k) * xyz_PiAdvSum(:,:,1) - xyz_PiAdvSum(:,:,k+1)
enddo
! $ \dot{\sigma} $ の上下境界値
! $ \dot{\sigma} $ on upper and lower boundary
!
xyr_SigDotN (:,:,0 ) = 0.
xyr_SigDotN (:,:,kmax) = 0.
xyr_SigDotNonG(:,:,0 ) = 0.
xyr_SigDotNonG(:,:,kmax) = 0.
! 温度の擾乱 (整数レベル), 仮温度, 仮温度の擾乱の計算
! Calculate temperature eddy (full level), virtual temperature,
! virtual temperature eddy
!
do k = 1, kmax
xyz_VirTemp(:,:,k) = &
& xyz_TempN(:,:,k) &
& * ( 1.0_DP + ((( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QVapN(:,:,k)) )
xyz_TempEdd(:,:,k) = xyz_TempN(:,:,k) - z_RefTemp(k)
xyz_VirTempEdd(:,:,k) = xyz_VirTemp(:,:,k) - z_RefTemp(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_TInpCoefA(k+1) * xyz_TempN(:,:,k+1) &
& + z_TInpCoefB(k ) * xyz_TempN(:,:,k ) &
& - r_RefTemp(k)
enddo
! 東西運動量移流項の計算
! Calculate advection of eastward momentum
!
xyz_UAdvN(:,:,1) = &
& ( xyz_VorN(:,:,1) + xy_Cori ) * xyz_VCosLatN(:,:,1) &
!
& - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) * xyr_SigDotN(:,:,1) &
& * ( xyz_UCosLatN(:,:,1) - xyz_UCosLatN(:,:,2) ) &
!
& - CpDry * z_TInpCoefK(1) * xyz_VirTempEdd(:,:,1) * xy_GradLambdaPiN
do k = 2, kmax-1
xyz_UAdvN(:,:,k) = &
& ( xyz_VorN(:,:,k) + xy_Cori ) * xyz_VCosLatN(:,:,k) &
!
& - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) &
& * ( &
& xyr_SigDotN(:,:,k-1) * ( xyz_UCosLatN(:,:,k-1) - xyz_UCosLatN(:,:,k) ) &
& + xyr_SigDotN(:,:,k) * ( xyz_UCosLatN(:,:,k) - xyz_UCosLatN(:,:,k+1) ) &
& ) &
!
& - CpDry * z_TInpCoefK(k) * xyz_VirTempEdd(:,:,k) * xy_GradLambdaPiN
end do
xyz_UAdvN(:,:,kmax) = &
& ( xyz_VorN(:,:,kmax) + xy_Cori ) * xyz_VCosLatN(:,:,kmax) &
!
& - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) * xyr_SigDotN(:,:,kmax-1) &
& * ( xyz_UCosLatN(:,:,kmax-1) - xyz_UCosLatN(:,:,kmax) ) &
!
& - CpDry * z_TInpCoefK(kmax) * xyz_VirTempEdd(:,:,kmax) * xy_GradLambdaPiN
! 南北運動量移流項の計算
! Calculate advection of northward momentum
!
xyz_VAdvN(:,:,1) = &
& - ( xyz_VorN(:,:,1) + xy_Cori ) * xyz_UCosLatN(:,:,1) &
!
& - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) * xyr_SigDotN(:,:,1)&
& * ( xyz_VCosLatN(:,:,1) - xyz_VCosLatN(:,:,2) ) &
!
& - CpDry * z_TInpCoefK(1) &
& * xyz_VirTempEdd(:,:,1) * xy_GradMuPiN
do k = 2, kmax-1
xyz_VAdvN(:,:,k) = &
& - ( xyz_VorN(:,:,k) + xy_Cori ) * xyz_UCosLatN(:,:,k) &
!
& - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) &
& * ( &
& xyr_SigDotN(:,:,k-1) * ( xyz_VCosLatN(:,:,k-1) - xyz_VCosLatN(:,:,k) ) &
& + xyr_SigDotN(:,:,k) * ( xyz_VCosLatN(:,:,k) - xyz_VCosLatN(:,:,k+1) ) &
& ) &
!
& - CpDry * z_TInpCoefK(k) * xyz_VirTempEdd(:,:,k) * xy_GradMuPiN
end do
xyz_VAdvN(:,:,kmax) = &
& - ( xyz_VorN(:,:,kmax) + xy_Cori ) * xyz_UCosLatN(:,:,kmax) &
!
& - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) * xyr_SigDotN(:,:,kmax-1) &
& * ( xyz_VCosLatN(:,:,kmax-1) - xyz_VCosLatN(:,:,kmax) ) &
!
& - CpDry * z_TInpCoefK(kmax) * xyz_VirTempEdd(:,:,kmax) * xy_GradMuPiN
! 運動エネルギー項 (仮温度補正含む) の計算
! Calculate kinematic energy term
! (including virtual temperature correction)
!
call HydroGrid( xyz_VirTemp - xyz_TempN, & ! (in)
& xyz_KinEngyN ) ! (out)
do k = 1, kmax
xyz_KinEngyN(:,:,k) = &
& xyz_KinEngyN(:,:,k) &
& + ( xyz_UCosLatN(:,:,k)**2 + xyz_VCosLatN(:,:,k)**2 ) &
& / ( 2.0_DP * ( 1.0_DP - xy_SinLat**2 ) )
end do
! 温度東西移流項, 温度南北移流項の計算
! Calculate eastward and northward advection of temperature
!
do k = 1, kmax
xyz_TempUAdvN(:,:,k) = xyz_UCosLatN(:,:,k) * xyz_TempEdd(:,:,k)
xyz_TempVAdvN(:,:,k) = xyz_VCosLatN(:,:,k) * xyz_TempEdd(:,:,k)
end do
! 温度の時間変化項 $ \hat{H} $ の計算
! Calculate temperature tendency term $ \hat{H} $
!
do k = 1, kmax-1
xyz_TempNonLinearN(:,:,k) = &
& xyz_TempEdd(:,:,k) * xyz_DivN(:,:,k) &
!
& - 1.0_DP / z_DelSigma(k) &
& * ( xyr_SigDotN(:,:,k-1) &
& * ( xyr_TempEdd(:,:,k-1) - xyz_TempEdd(:,:,k) ) &
& + xyr_SigDotN(:,:,k) &
& * ( xyz_TempEdd(:,:,k) - xyr_TempEdd(:,:,k) ) ) &
!
& - 1.0_DP / z_DelSigma(k) &
& * ( xyr_SigDotNonG(:,:,k-1) &
& * ( r_RefTemp(k-1) - z_RefTemp(k) ) &
& + xyr_SigDotNonG(:,:,k) &
& * ( z_RefTemp(k) - r_RefTemp(k) ) ) &
!
& + z_TInpCoefK(k) * xyz_VirTemp(:,:,k) * xyz_PiAdv(:,:,k) &
!
& - z_HydroAlpha(k) / z_DelSigma(k) &
& * ( xyz_VirTemp(:,:,k) * xyz_PiAdvSum(:,:,k) &
& + xyz_VirTempEdd(:,:,k) * xyz_DivSum(:,:,k) ) &
!
& - z_HydroBeta(k) / z_DelSigma(k) &
& * ( xyz_VirTemp(:,:,k) * xyz_PiAdvSum(:,:,k+1) &
& + xyz_VirTempEdd(:,:,k) * xyz_DivSum(:,:,k+1) )
enddo
xyz_TempNonLinearN(:,:,kmax) = &
& xyz_TempEdd(:,:,kmax) * xyz_DivN(:,:,kmax) &
!
& - 1.0_DP / z_DelSigma(kmax) &
& * ( xyr_SigDotN(:,:,kmax-1) &
& * ( xyr_TempEdd(:,:,kmax-1) - xyz_TempEdd(:,:,kmax) ) &
& + xyr_SigDotN(:,:,kmax) &
& * ( xyz_TempEdd(:,:,kmax) - xyr_TempEdd(:,:,kmax) ) ) &
!
& - 1.0_DP / z_DelSigma(kmax) &
& * ( xyr_SigDotNonG(:,:,kmax-1) &
& * ( r_RefTemp(kmax-1) - z_RefTemp(kmax) ) &
& + xyr_SigDotNonG(:,:,kmax) &
& * ( z_RefTemp(kmax) - r_RefTemp(kmax) ) ) &
!
& + z_TInpCoefK(kmax) * xyz_VirTemp(:,:,kmax) * xyz_PiAdv(:,:,kmax) &
!
& - z_HydroAlpha(kmax) / z_DelSigma(kmax) &
& * ( xyz_VirTemp(:,:,kmax) * xyz_PiAdvSum(:,:,kmax) &
& + xyz_VirTempEdd(:,:,kmax) * xyz_DivSum(:,:,kmax) )
end subroutine NonLinearOnGrid
!--------------------------------------------------------------------------------------
subroutine NonLinearOnGridQMix( &
& xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, & ! (in)
& xyzf_QMixN, & ! (in)
!
& xyzf_QMixNonLinearN, & ! (out)
& xyzf_QMixUAdvN, xyzf_QMixVAdvN & ! (out)
& )
!
! 非線形項 (非重力波項) を格子点上で計算します.
!
! Non-linear terms (non gravitational terms) are calculated on
! grid points
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! 格子点設定
! Grid points settings
!
use gridset, only: imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax
! 成分の数
! Number of composition
implicit none
real(DP), intent(in):: xyz_UCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ U (t) = u (t) \cos \varphi $ .
real(DP), intent(in):: xyz_VCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ V (t) = v (t) \cos \varphi $ .
real(DP), intent(in):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
! $ D (t) $ . 発散. Divergence
real(DP), intent(in ):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} (t) $ .
! 鉛直流. Vertical flow
real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t) $ . 比湿. Specific humidity
real(DP), intent(out):: xyzf_QMixNonLinearN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ R (t) $ . 比湿時間変化項.
! Specific humidity tendency
real(DP), intent(out):: xyzf_QMixUAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Uq (t) $ . 比湿東西移流項.
! Eastward advection of specific humidity
real(DP), intent(out):: xyzf_QMixVAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Vq (t) $ . 比湿南北移流項.
! Northward advection of specific humidity
!-----------------------------------
! 作業変数
! 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
!
! 比湿東西移流項, 比湿南北移流項の計算
! Calculate eastward and northward advection of specific humidity
!
do n = 1, ncmax
do k = 1, kmax
xyzf_QMixUAdvN(:,:,k,n) = xyz_UCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
xyzf_QMixVAdvN(:,:,k,n) = xyz_VCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
end do
end do
! 比湿時間変化項 $ R $ の計算
! Calculate specific humidity tendency $ R $
!
do n = 1, ncmax
xyzf_QMixNonLinearN(:,:,1,n) = &
& xyzf_QMixN(:,:,1,n) * xyz_DivN(:,:,1) &
& - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) &
& * xyr_SigDotN(:,:,1) &
& * ( xyzf_QMixN(:,:,1,n) - xyzf_QMixN(:,:,2,n) )
do k = 2, kmax - 1
xyzf_QMixNonLinearN(:,:,k,n) = &
& xyzf_QMixN(:,:,k,n) * xyz_DivN(:,:,k) &
& - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) &
& * ( &
& xyr_SigDotN(:,:,k-1) &
& * ( xyzf_QMixN(:,:,k-1,n) - xyzf_QMixN(:,:,k ,n) ) &
& + xyr_SigDotN(:,:,k) &
& * ( xyzf_QMixN(:,:,k ,n) - xyzf_QMixN(:,:,k+1,n) ) )
end do
xyzf_QMixNonLinearN(:,:,kmax,n) = &
& xyzf_QMixN(:,:,kmax,n) * xyz_DivN(:,:,kmax) &
& - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) &
& * xyr_SigDotN(:,:,kmax-1) &
& * ( xyzf_QMixN(:,:,kmax-1,n) - xyzf_QMixN(:,:,kmax,n) )
end do
end subroutine NonLinearOnGridQMix
!--------------------------------------------------------------------------------------
subroutine HydroGrid( xyz_Temp, & ! (in)
& xyz_Phi & ! (out)
& )
!
! 格子点データである温度 $ T $ から, 静水圧の式を用いて
! 格子点データのジオポテンシャル高度 $ \Phi $ を求めます.
!
! モジュール引用 ; USE statements
!
use constants, only: &
& CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 格子点設定
! Grid points settings
!
use gridset, only: imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 宣言文 ; 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( &
& w_SurfGeoPot, & ! (in)
& wz_DVorDtNG, wz_DDivDtNG, wz_DTempDtNG, w_DPiDtNG, & ! (in)
& wz_DivN, wz_TempN, w_PiN, & ! (in)
& wz_VorB, wz_DivB, wz_TempB, w_PiB, & ! (in)
& wz_VorA, wz_DivA, wz_TempA, w_PiA & ! (out)
& )
!
! 時間積分を行い,
! 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から
! 時刻 $ t+\Delta t $ の物理量を計算します.
!
! 時間積分法にはリープフロッグスキームを用いています.
! ただし拡散項には後方差分を用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! With time integration, physical values at $ t+\Delta t $ is calculated
! from tendency at $ t $ and physical values at $ t-\Delta t $ .
!
! Leap-frog scheme is used as time integration scheme.
! And backward scheme is used for diffusion terms.
! As default setting, semi-implicit scheme is applied to gravitational terms
! in order to enlarge the value 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, &
! $ a $ [m].
! 惑星半径.
! Radius of planet
& CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! LU 分解法により連立 1 次方程式を解くための関数 (spml 同梱モジュール)
! Functions to solve linear simultaneous equation by LU decomposition
! (a module included in spml)
!
use lumatrix, only: LUSolve
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax
! 成分の数
! Number of composition
implicit none
real(DP), intent(in):: w_SurfGeoPot (lmax)
! $ \Phi_s $ . 地表ジオポテンシャル.
! Surface geo-potential
real(DP), intent(in):: wz_DVorDtNG (lmax, 1:kmax)
! $ \left( \DD{\zeta}{t} (t) \right)^{NG}$ . 渦度変化の非重力波成分 (スペクトル).
! Non-gravity wave component of vorticity tendency (spectral)
real(DP), intent(in):: wz_DDivDtNG (lmax, 1:kmax)
! $ \DD{D}{t} (t) $ . 発散変化の非重力波成分 (スペクトル).
! Divergence tendency (spectral)
real(DP), intent(in):: wz_DTempDtNG(lmax, 1:kmax)
! $ \left( \DD{T}{t} (t) \right)^{NG}$ . 温度変化の非重力波成分 (スペクトル).
! Temperature tendency (spectral)
real(DP), intent(in):: w_DPiDtNG(lmax)
! $ \left( \DD{p_s}{t} (t) \right) $ . 地表面気圧変化の非重力波項 (スペクトル).
! Non-gravity wave component of surface pressure tendency (spectral)
real(DP), intent(in):: wz_DivN (lmax, 1:kmax)
! $ D (t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(in):: wz_TempN (lmax, 1:kmax)
! $ T (t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(in):: w_PiN (lmax)
! $ \pi = \ln p_s (t) $ . 地表面気圧 (スペクトル).
real(DP), intent(in):: wz_VorB (lmax, 1:kmax)
! $ \zeta (t-\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(in):: wz_DivB (lmax, 1:kmax)
! $ D (t-\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(in):: wz_TempB (lmax, 1:kmax)
! $ T (t-\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(in):: w_PiB (lmax)
! $ \pi = \ln p_s (t-\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP), intent(out):: wz_VorA (lmax, 1:kmax)
! $ \zeta (t+\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(out):: wz_DivA (lmax, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(out):: wz_TempA (lmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(out):: w_PiA (lmax)
! $ \pi = \ln p_s (t+\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
! 作業変数
! Work variables
!
real(DP):: wz_HDiv (lmax, 1:kmax)
real(DP):: w_CtDiv (lmax)
real(DP):: wz_WT (lmax, 1:kmax)
real(DP):: wz_siTemp (lmax, 1:kmax)
! 温度 (セミインプリシット法のための作業変数).
! Temperature (work variable for semi-implicit scheme)
real(DP):: w_siPi (lmax)
! $ \pi $ (セミインプリシット法のための作業変数).
! $ \pi $ (work variable for semi-implicit scheme)
real(DP):: wz_siPhi (lmax, 1:kmax)
! $ \Phi = \underline{W} ( 1 - 2 \Delta t \underline{D_H} ) \overline{ \Dvect{T} }^{t}$ .
! (セミインプリシット法のための作業変数).
! (Work variable for semi-implicit scheme)
real(DP):: wz_siVectF (lmax, 1:kmax)
! $ \Dvect{f} $ .
! 発散項に関するセミインプリシット方程式の右辺.
! Right-hand side of a semi-implicit equation of a divergence term.
real(DP):: wz_siDivAvrTime (lmax, 1:kmax)
! $ \overline{\Dvect{D}}^{t} $ .
! 時間平均の $ \Dvect{D} $ (セミインプリシット法のための作業変数).
! Time average $ \Dvect{D} $ (a work variable for semi-implicit scheme)
integer:: k, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
! 時間積分. 拡散は後方差分
! Time integration. Backward difference is applied to diffusion
!
! 渦度 ; Vorticity
!
wz_VorA = &
& ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM ) ) &
& * ( wz_VorB + 2.0_DP * DelTime * wz_DVorDtNG )
! 発散 ; Divergence
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeSemiImplicit )
! セミインプリシット法で用いる行列の計算に使われる項の計算
! Calculate terms used in making a matrix for semi-implicit method
!
wz_siTemp = &
& ( 1.0_DP - DelTime * wz_DisCoefH ) * wz_TempB &
& + DelTime * wz_DTempDtNG
! NOTE:
! The matrix wz_siMtxDi = ( 1 - 2 \Delta t D_H )^{-1} is diagonal matrix.
!
wz_siPhi (:,1) = CpDry * z_HydroAlpha(1) * wz_siMtxDi(:,1) * wz_siTemp(:,1)
do k = 2, kmax
wz_siPhi (:,k) = wz_siPhi(:,k-1) &
& + CpDry * z_HydroAlpha(k ) * wz_siMtxDi(:,k ) * wz_siTemp(:,k ) &
& + CpDry * z_HydroBeta (k-1) * wz_siMtxDi(:,k-1) * wz_siTemp(:,k-1)
end do
w_siPi = w_PiB + DelTime * w_DPiDtNG
! 発散方程式の右辺の計算
! Calculate right side of divergence equation
!
do k = 1, kmax
wz_siVectF(:,k) = &
& ( 1.0_DP - DelTime * wz_DisCoefM (:,k) ) * wz_DivB(:,k) &
!
& + DelTime * wz_DDivDtNG(:,k) &
!
& - DelTime * w_LaplaEigVal(:) &
& * ( w_SurfGeoPot + wz_siPhi(:,k) + z_siMtxG(k) * w_siPi )
end do
! 時間平均の $ \Dvect{D} $ を LU 行列で解く
! Solve time average $ \Dvect{D} $ with LU matrix
!
wz_siDivAvrTime = LUSolve( wzz_siMtxLU, wz_siMtxPiv, wz_siVectF )
wz_DivA = 2. * wz_siDivAvrTime - wz_DivB
case ( IDTimeIntegSchemeExplicit )
wz_WT = 0.0_DP
do k = 1, kmax
do kk = 1, kmax
wz_WT(:,k) = wz_WT(:,k) + zz_siMtxW(k,kk) * wz_TempN(:,kk)
end do
end do
do k = 1, kmax
wz_DivA(:,k) = &
& ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM(:,k) ) ) &
& * ( wz_DivB(:,k) &
& + 2.0_DP * DelTime &
& * ( wz_DDivDtNG(:,k) &
& - w_LaplaEigVal(:) &
& * ( w_SurfGeoPot &
& + wz_WT(:,k) &
& + z_siMtxG(k) * w_PiN &
& ) &
& ) &
& )
end do
end select
! 温度 ; Temperature
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeSemiImplicit )
wz_HDiv = 0.
do k = 1, kmax
do kk = 1, kmax
wz_HDiv(:,k) = wz_HDiv(:,k) + zz_siMtxH(k,kk) * wz_siDivAvrTime(:,kk)
end do
end do
case ( IDTimeIntegSchemeExplicit )
wz_HDiv = 0.
do k = 1, kmax
do kk = 1, kmax
wz_HDiv(:,k) = wz_HDiv(:,k) + zz_siMtxH(k,kk) * wz_DivN(:,kk)
end do
end do
end select
wz_TempA = &
& ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH ) ) &
& * ( wz_TempB + 2.0_DP * DelTime * ( wz_DTempDtNG - wz_HDiv ) )
! 地表面気圧 ; Surface pressure
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeSemiImplicit )
w_CtDiv = 0.0_DP
do k = 1, kmax
w_CtDiv = w_CtDiv + z_siMtxC(k) * wz_siDivAvrTime(:,k)
end do
case ( IDTimeIntegSchemeExplicit )
w_CtDiv = 0.0_DP
do k = 1, kmax
w_CtDiv = w_CtDiv + z_siMtxC(k) * wz_DivN(:,k)
end do
end select
w_PiA = w_PiB + 2.0_DP * DelTime * ( w_DPiDtNG - w_CtDiv )
end subroutine TimeIntegration
!--------------------------------------------------------------------------------------
subroutine TimeIntegrationQMix( &
& wzf_DQMixDtN, & ! (in)
& wzf_QMixB, & ! (in)
& wzf_QMixA & ! (out)
& )
!
! 時間積分を行い,
! 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から
! 時刻 $ t+\Delta t $ の物理量を計算します.
!
! 時間積分法にはリープフロッグスキームを用いています.
! ただし拡散項には後方差分を用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! With time integration, physical values at $ t+\Delta t $ is calculated
! from tendency at $ t $ and physical values at $ t-\Delta t $ .
!
! Leap-frog scheme is used as time integration scheme.
! And backward scheme is used for diffusion terms.
! As default setting, semi-implicit scheme is applied to gravitational terms
! in order to enlarge the value of $ \Delta t $ .
! Explicit scheme can be applied to gravitational terms by changing
! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
!
! モジュール引用 ; USE statements
!
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax
! 成分の数
! Number of composition
implicit none
real(DP), intent(in):: wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax)
! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP), intent(in):: wzf_QMixB(lmax, 1:kmax, 1:ncmax)
! $ q (t-\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP), intent(out):: wzf_QMixA(lmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
! 作業変数
! 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
!
! 時間積分. 拡散は後方差分
! Time integration. Backward difference is applied to diffusion
!
! 比湿 ; Specific humidity
!
do n = 1, ncmax
do k = 1, kmax
wzf_QMixA(:,k,n) = &
& ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * w_DisCoefQ ) ) &
& * ( wzf_QMixB(:,k,n) + 2.0_DP * DelTime * wzf_DQMixDtN(:,k,n) )
end do
end do
end subroutine TimeIntegrationQMix
!--------------------------------------------------------------------------------------
subroutine DivergenceDamping( &
& wz_DivA & ! (inout)
& )
!
! 発散の減衰項の付加 (場がつり合っていない計算初期の擾乱の減衰を想定)
!
! Addition of divergence damping term
!
! モジュール引用 ; USE statements
!
! 時刻管理
! Time control
!
use timeset, only: &
& DelTime, & ! $ \Delta t $ [s]
& TimeN ! ステップ $ t $ の時刻. Time of step $ t $.
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
implicit none
real(DP), intent(inout):: wz_DivA (lmax, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
! 作業変数
! Work variables
!
real(DP) :: DivDampCoef
! 実行文 ; Executable statement
!
if ( FlagDivDamp ) then
DivDampCoef = 1.0_DP / DelTime * ( DivDampPeriod - ( TimeN - DivDampTime0 ) ) / DivDampPeriod
if ( DivDampCoef < 0.0_DP ) DivDampCoef = 0.0_DP
wz_DivA = wz_DivA / ( 1.0_DP + 2.0_DP * DelTime * DivDampCoef )
end if
end subroutine DivergenceDamping
!--------------------------------------------------------------------------------------
subroutine OutputDiagnosedVariables( &
& xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, & ! (in)
& xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, & ! (in)
& xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA, & ! (in)
& xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, & ! (in)
& xyr_SigDot, xy_DPiDt, & ! (in)
& xyz_PiAdv, & ! (in)
& xyz_OMG & ! (out)
& )
!
! 診断量の出力を行います.
!
! Diagnostic variables are output.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: &
& RPlanet, &
! $ a $ [m].
! 惑星半径.
! Radius of planet
& Grav, &
! $ g $ [m s-2].
! 重力加速度.
! Gravitational acceleration
& CpDry, &
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
& LatentHeat
! $ L $ [J kg-1] .
! 凝結の潜熱.
! Latent heat of condensation
! 時刻管理
! Time control
!
use timeset, only: &
& DelTime, & ! $ \Delta t $ [s]
& TimeN ! ステップ $ t $ の時刻. Time of step $ t $.
! 格子点設定
! Grid points settings
!
use gridset, only: &
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax, &
! 成分の数
! Number of composition
& a_QMixName, &
& IndexH2OVap
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& z_Sigma, & ! $ \sigma $ レベル (整数).
! Full $ \sigma $ level
& r_Sigma, & ! $ \sigma $ レベル (半整数).
! Half $ \sigma $ level
& z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
!
! 積分と平均の操作
! Operation for integral and average
!
use intavr_operate, only: &
IntLonLatSig_xyz, AvrLonLatSig_xyz, &
IntLonLat_xy, AvrLonLat_xy
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef AXISYMMETRY
use wa_zonal_module, only: &
& w_xy, &
& xy_GradLon_w, xy_GradLat_w, &
& wa_DivLambda_xya, &
& wa_DivMu_xya, &
& xya_wa
#elif LIB_MPI
! modify for bug of Intel fortran (2025/09/11 S.Takehiro)
! use ua_mpi_module, only: &
! & w_xy => u_pv, &
! & xy_GradLon_w => pv_GradLon_u, xy_GradLat_w => pv_GradLat_u, &
! & wa_DivLambda_xya => ua_DivLambda_pva, &
! & wa_DivMu_xya => ua_DivMu_pva, &
! & xya_wa => pva_ua
use ua_mpi_module_base, only: &
& w_xy => u_pv, &
& xya_wa => pva_ua
use ua_mpi_module_deriv, only: &
& xy_GradLon_w => pv_GradLon_u, xy_GradLat_w => pv_GradLat_u, &
& wa_DivLambda_xya => ua_DivLambda_pva, &
& wa_DivMu_xya => ua_DivMu_pva
#else
use wa_module, only: &
& w_xy, &
& xy_GradLon_w, xy_GradLat_w, &
& wa_DivLambda_xya, &
& wa_DivMu_xya, &
& xya_wa
#endif
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 宣言文 ; 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):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ 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):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ 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_UA (0:imax-1, 1:jmax, 1:kmax)
! $ u (t+\Delta t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VA (0:imax-1, 1:jmax, 1:kmax)
! $ v (t+\Delta t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度. Temperature
real(DP), intent(in):: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsA (0:imax-1, 1:jmax)
! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure
real(DP), intent(in):: xyz_DUDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{u}{t}\right)^{phy} $ .
! 外力項 (物理過程) による東西風速変化.
! Eastward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DVDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{v}{t}\right)^{phy} $ .
! 外力項 (物理過程) による南北風速変化.
! Northward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DTempDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{T}{t}\right)^{phy} $ .
! 外力項 (物理過程) による温度変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \left(\DP{q}{t}\right)^{phy} $ .
! 外力項 (物理過程) による比湿変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xyr_SigDot (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流. Vertical flow
real(DP), intent(in):: xy_DPiDt (0:imax-1, 1:jmax)
! $ Z $ . 地表面気圧時間変化項.
! Surface pressure tendency
real(DP), intent(in):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{v} \cdot \nabla \pi $ .
! $ \pi $ の移流. Advection of $ \pi $
real(DP), intent(out):: xyz_OMG (0:imax-1, 1:jmax, 1:kmax)
!
! OMEGA, DP/Dt
! 作業変数
! Work variables
!
real(DP):: xyz_DUDtDyn (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_DVDtDyn (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_DTempDtDyn (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyzf_DQMixDtDyn(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP):: xy_DPsDtDyn (0:imax-1, 1:jmax)
real(DP):: xyz_UCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ U = u \cos \varphi $ .
real(DP):: xyz_VCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ V = v \cos \varphi $ .
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):: xyr_PiAdvDivSum(0:imax-1, 1:jmax, 0:kmax)
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
integer:: n
! 実行文 ; Executable statement
!
! Calculate tendencies
!
xyz_DUDtDyn = ( xyz_UA - xyz_UB ) / ( 2.0_DP * DelTime ) &
& - xyz_DUDtPhy
xyz_DVDtDyn = ( xyz_VA - xyz_VB ) / ( 2.0_DP * DelTime ) &
& - xyz_DVDtPhy
xyz_DTempDtDyn = ( xyz_TempA - xyz_TempB ) / ( 2.0_DP * DelTime ) &
& - xyz_DTempDtPhy
xyzf_DQMixDtDyn = ( xyzf_QMixA - xyzf_QMixB ) / ( 2.0_DP * DelTime ) &
& - xyzf_DQMixDtPhy
xy_DPsDtDyn = ( xy_PsA - xy_PsB ) / ( 2.0_DP * DelTime )
call HistoryAutoPut( TimeN, 'DUDtDyn' , xyz_DUDtDyn )
call HistoryAutoPut( TimeN, 'DVDtDyn' , xyz_DVDtDyn )
call HistoryAutoPut( TimeN, 'DTempDtDyn', xyz_DTempDtDyn )
do n = 1, ncmax
call HistoryAutoPut( TimeN, 'D'//trim(a_QMixName(n))//'DtDyn', xyzf_DQMixDtDyn(:,:,:,n) )
end do
call HistoryAutoPut( TimeN, 'DPsDtDyn' , xy_DPsDtDyn )
! 鉛直流と地表面気圧時間変化項の出力
! Output vertical flow and surface pressure tendency
!
call HistoryAutoPut( TimeN, 'SigDot', xyr_SigDot )
call HistoryAutoPut( TimeN, 'DPiDt', xy_DPiDt )
! 風速から渦度発散の計算
! Calculate vorticity and divergence from wind velocity
!
do k = 1, kmax
xyz_UCosLat(:,:,k) = xyz_UN(:,:,k) * xy_CosLat
xyz_VCosLat(:,:,k) = xyz_VN(:,:,k) * xy_CosLat
end do
xyz_Vor = xya_wa( wa_DivLambda_xya( xyz_VCosLat ) &
& - wa_DivMu_xya( xyz_UCosLat ) ) / RPlanet
xyz_Div = xya_wa( wa_DivLambda_xya( xyz_UCosLat ) &
& + wa_DivMu_xya( xyz_VCosLat ) ) / RPlanet
call HistoryAutoPut( TimeN, 'Vor', xyz_Vor )
call HistoryAutoPut( TimeN, 'Div', xyz_Div )
! Calculation of Omega (DPressDt)
! It should be recognized that the value of xy_Ps here may have been modified
! by mass fixer.
!
! Integration from top of the model to k's layer upper interface
xyr_PiAdvDivSum(:,:,kmax) = 0.0_DP
do k = kmax-1, 0, -1
xyr_PiAdvDivSum(:,:,k) = xyr_PiAdvDivSum(:,:,k+1) &
& + ( xyz_PiAdv(:,:,k+1) + xyz_Div(:,:,k+1) ) * z_DelSigma(k+1)
end do
do k = 1, kmax
xyz_OMG(:,:,k) = &
& xy_PsN * ( &
& z_Sigma(k) * xyz_PiAdv(:,:,k) &
& - xyr_PiAdvDivSum(:,:,k) &
& - ( xyz_PiAdv(:,:,k) + xyz_Div(:,:,k) ) * ( z_Sigma(k) - r_Sigma(k) ) &
& )
end do
call HistoryAutoPut( TimeN, 'OMG', xyz_OMG )
! 質量の計算
! Calculate mass
!
xy_Mass = xy_PsN / Grav
! エネルギー, エンストロフィーの計算
! Calculate energy and enstrophy
!
call HydroGrid( xyz_TempN, & ! (in)
& xyz_PotEngy & ! (out)
& )
do k = 1, kmax
xyz_KinEngy(:,:,k) = ( xyz_UN(:,:,k) ** 2 + xyz_VN(:,:,k) ** 2 ) / 2.0_DP &
& * xy_Mass
xyz_IntEngy(:,:,k) = CpDry * xyz_TempN(:,:,k) &
& * xy_Mass
xyz_PotEngy(:,:,k) = xyz_PotEngy(:,:,k) &
& * xy_Mass
xyz_LatEngy(:,:,k) = LatentHeat * xyzf_QMixN(:,:,k,IndexH2OVap) &
& * 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 )
! 空間平均された診断量
! Spatially averaged diagnostic variables
!
call HistoryAutoPut( TimeN, 'TotalMass', IntLonLat_xy(xy_Mass) )
call HistoryAutoPut( TimeN, 'TotalKinEngy', IntLonLatSig_xyz(xyz_KinEngy) )
call HistoryAutoPut( TimeN, 'TotalIntEngy', IntLonLatSig_xyz(xyz_IntEngy) )
call HistoryAutoPut( TimeN, 'TotalPotEngy', IntLonLatSig_xyz(xyz_PotEngy) )
call HistoryAutoPut( TimeN, 'TotalLatEngy', IntLonLatSig_xyz(xyz_LatEngy) )
call HistoryAutoPut( TimeN, 'TotalTotEngy', IntLonLatSig_xyz(xyz_TotEngy) )
call HistoryAutoPut( TimeN, 'TotalEnstro', IntLonLatSig_xyz(xyz_Enstro) )
call HistoryAutoPut( TimeN, 'MeanMass', AvrLonLat_xy(xy_Mass) )
call HistoryAutoPut( TimeN, 'MeanKinEngy', AvrLonLatSig_xyz(xyz_KinEngy) )
call HistoryAutoPut( TimeN, 'MeanIntEngy', AvrLonLatSig_xyz(xyz_IntEngy) )
call HistoryAutoPut( TimeN, 'MeanPotEngy', AvrLonLatSig_xyz(xyz_PotEngy) )
call HistoryAutoPut( TimeN, 'MeanLatEngy', AvrLonLatSig_xyz(xyz_LatEngy) )
call HistoryAutoPut( TimeN, 'MeanTotEngy', AvrLonLatSig_xyz(xyz_TotEngy) )
call HistoryAutoPut( TimeN, 'MeanEnstro', AvrLonLatSig_xyz(xyz_Enstro) )
end subroutine OutputDiagnosedVariables
!--------------------------------------------------------------------------------------
subroutine DynamicsHSplVAS83Init
!
! 計算に必要なパラメタの設定や NAMELIST#dynamics_hspl_vas83_nml
! の読み込みを行います.
!
! Configure parameters for calculation,
! and load "NAMELIST#dynamics_hspl_vas83_nml"
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: &
& RPlanet, &
! $ a $ [m].
! 惑星半径.
! Radius of planet
& Omega, &
! $ \Omega $ [s-1].
! 回転角速度.
! Angular velocity
& GasRDry, &
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
& CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 格子点設定
! Grid points settings
!
use gridset, only: nmax, & ! 最大全波数.
! Maximum truncated wavenumber
& lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax, &
! 成分の数
! Number of composition
& a_QMixName
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& z_Sigma, & ! $ \sigma $ レベル (整数).
! Full $ \sigma $ level
& r_Sigma, & ! $ \sigma $ レベル (半整数).
! Half $ \sigma $ level
& z_DelSigma, & ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
& AxNameX, AxNameY, AxNameZ, AxNameR, AxNameT
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef AXISYMMETRY
use wa_zonal_module, only: &
& xy_Lat, &
& w_xy, &
& rn, &
& nm_l
#elif LIB_MPI
use ua_mpi_module, only: &
& xy_Lat => pv_Lat, &
& w_xy => u_pv, &
& rn => rnu, &
& nm_l => nm_lu
#else
use wa_module, only: &
& xy_Lat, &
& w_xy, &
& rn, &
& nm_l
#endif
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg
! メッセージ制御
! Message control
!
use mpi_messagecntl, only : DoesOutputMPIMessage
! ヒストリデータ出力
! 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, & ! 標準出力の装置番号. Unit number of standard output
& STRING ! 文字列. Strings.
! 日付および時刻の取り扱い
! Date and time handler
!
use dc_calendar, only: DCCalConvertByUnit
! 組み込み関数 PRESENT の拡張版関数
! Extended functions of intrinsic function "PRESENT"
!
use dc_present, only: present_and_not_empty
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 時刻管理
! Time control
!
use timeset, only: &
& TimeN ! ステップ $ t $ の時刻. Time of step $ t $.
!
! 水平スペクトル解析
! horizontal spectral analysis
!
use dynamics_diag_spectrum, only : DynamicsDiagSpAnalysisInit
! 質量の補正
! Mass fixer
!
use mass_fixer, only : MassFixerInit
! セミラグランジュ法による物質移流計算
! Semi-Lagrangian method for tracer transport
use sltt, only: SLTTInit
!!$ ! 物質移流 (有限体積法)
!!$ ! Tracer Transport (finite volume method)
!!$ !
!!$ use tt_fv, only: TTFVInit
!
! Utility module for advection test
!
use adv_test, only: AdvTestInit
! 宣言文 ; Declaration statements
!
implicit none
! 基準温度の設定のための作業変数
! Work variable for reference temperature
!
character(TOKEN) :: TimeIntegScheme
! 時間積分法.
! 以下の方法を選択可能.
!
! Time integration scheme.
! Available schemes are as follows.
!
! * "Semi-implicit"
! * "Explicit"
real(DP):: RefTemp
! $ \overline{T} $ . 基準温度.
! Reference temperature
! 水平拡散, スポンジ層のための作業変数
! Work variable for horizontal diffusion and sponge layer
!
real(DP) :: w_HDifCoefM (lmax)
! $ {\cal D}_M = -K_{HD} [(-1)^{N_D/2}\nabla^{N_D}- (2/a^2)^{N_D/2}] $ .
! 運動量水平拡散係数.
! Coefficient of horizontal momentum diffusion
real(DP) :: w_HDifCoefH (lmax)
! $ {\cal D}_H = -(-1)^{N_D/2}K_{HD}\nabla^{N_D} $ .
! 熱, 水水平拡散係数.
! Coefficient of horizontal thermal and water diffusion
real(DP) :: wz_SpongeLayerCoefM(lmax, 1:kmax)
! スポンジ層における運動量の減衰係数
! Damping coefficient for momentum in a sponge layer
real(DP) :: wz_SpongeLayerCoefH(lmax, 1:kmax)
! スポンジ層における熱の減衰係数
! Damping coefficient for temperature zonal perturbation
! in a sponge layer
integer:: HDOrder ! 超粘性の次数. Order of hyper-viscosity
real(DP):: VisCoef ! 超粘性係数. Hyper-viscosity coefficient
real(DP):: HDEFoldTimeValue ! 最大波数に対する e-folding time.
! 負の値を与えると, 水平拡散係数をゼロにします.
!
! E-folding time for maximum wavenumber.
! If negative value is given,
! coefficients of horizontal diffusion become zero.
character(TOKEN):: HDEFoldTimeUnit
! 最大波数に対する e-folding time の単位.
! Unit of e-folding time for maximum wavenumber
real(DP):: HDEFoldTime ! 最大波数に対する e-folding time [単位: 秒].
! E-folding time for maximum wavenumber [Unit: sec]
logical :: FlagSpongeLayer
logical :: FlagSpongeLayerforZonalMean
logical :: FlagSpongeLayerforHeat
real(DP) :: SLEFoldTimeValue
! スポンジ層の最上層における減衰時定数
! Damping time constant at the top model layer
! in a sponge layer
character(TOKEN) :: SLEFoldTimeUnit
! SLEFoldTimeValue の単位
! Unit of SLEFoldTimeValue
real(DP) :: SLEFoldTime
! スポンジ層の最上層における減衰時定数 (秒)
! Damping time constant at the top model layer
! in a sponge layer (sec)
integer :: SLOrder
! スポンジ層の減衰係数の sigma 依存性のオーダ
! Sigma dependence of damping coefficient
! in a sponge layer
integer :: SLNumLayer
! スポンジ層が適応されるモデル上端からの層の数
! Number of layers which the sponge layer is applied.
integer :: a_DegOrd(2)
real(DP) :: DivDampPeriodValue
! Period for divergence damping application
character(TOKEN) :: DivDampPeriodUnit
! Unit of DivDampPeriodValue
!!$ character(TOKEN):: TracerTransportMethod
! NonLinearOnGrid 等で使用する係数の設定のための作業変数
! Work variable for coefficients for "NonLinearOnGrid", etc.
!
real(DP):: Kappa ! $ \kappa = R/C_p $ .
! 気体定数の定圧比熱に対する比.
! Ratio of gas constant to specific heat
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n
integer:: l ! 波数方向に回る DO ループ用作業変数
! Work variables for DO loop in wavenumber
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, &
& HDOrder, &
& HDEFoldTimeValue, HDEFoldTimeUnit, &
& FlagSpongeLayer, FlagSpongeLayerforZonalMean, FlagSpongeLayerforHeat, &
& SLEFoldTimeValue, SLEFoldTimeUnit, SLOrder, SLNumLayer, &
& RefTemp, &
& FlagDivDamp, DivDampPeriodValue, DivDampPeriodUnit, &
& DivDampTime0, &
& FlagTotMassFix, &
& FlagMassHorDifCor, &
!!$ & TracerTransportMethod, &
& FlagSLTT, &
& FlagAdvTest
!
! デフォルト値については初期化手続 "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
! デフォルト値の設定
! Default values settings
!
TimeIntegScheme = 'Semi-implicit'
HDOrder = 8
HDEFoldTimeValue = 8640.0_DP
HDEFoldTimeUnit = 'sec'
FlagSpongeLayer = .false.
FlagSpongeLayerforZonalMean = .false.
FlagSpongeLayerforHeat = .false.
SLEFoldTimeValue = 86400.0_DP
SLEFoldTimeUnit = 'sec'
SLOrder = 1
SLNumLayer = kmax
RefTemp = 300.
FlagDivDamp = .false.
DivDampPeriodValue = 2.0_DP
DivDampPeriodUnit = 'day'
DivDampTime0 = 0.0_DP
FlagTotMassFix = .true.
FlagMassHorDifCor = .false.
!!$ TracerTransportMethod = "spectrum"
FlagSLTT = .false.
FlagAdvTest = .false.
! NAMELIST の読み込み
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, & ! (out)
& namelist_filename, mode = 'r' ) ! (in)
rewind( unit_nml )
read( unit_nml, & ! (in)
& nml = dynamics_hspl_vas83_nml, & ! (out)
& iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = dynamics_hspl_vas83_nml )
end if
! 時間積分法のチェック
! Check time integration scheme
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
IDTimeIntegScheme = IDTimeIntegSchemeSemiImplicit
case ('explicit')
IDTimeIntegScheme = IDTimeIntegSchemeExplicit
case default
call MessageNotify( 'E', module_name, &
& '"TimeIntegScheme" must be "Semi-Implicit" or "Explicit".' )
end select
!
! Check tracer transport method
!
!!$ select case ( LChar( trim( TracerTransportMethod ) ) )
!!$ case ('spectrum')
!!$ IDTTMethod = IDTTMethodSp
!!$ case ('semi-lagrange')
!!$ IDTTMethod = IDTTMethodSL
!!$ case ('finite volume')
!!$ IDTTMethod = IDTTMethodFV
!!$ case default
!!$ call MessageNotify( 'E', module_name, &
!!$ & '"TracerTransportMethod" must be "spectrum", "semi-lagrange", or "finite volume".' )
!!$ end select
! SemiImplMatrix サブルーチン用に $ \Delta t $ の保存値に初期値を設定.
! Configure initial value to saved value of $ \Delta t $ for a subroutine "SemiImplMatrix"
!
DelTimeSave = - 1.0_DP
! NonLinearOnGrid 等で使用する係数の設定
! Configure coefficients for "NonLinearOnGrid", etc.
!
! $ \sin \varphi $ と $ \cos \varphi $ の計算
! Calculate $ \sin \varphi $ and $ \cos \varphi $
!
allocate( xy_SinLat (0:imax-1, 1:jmax) )
allocate( xy_CosLat (0:imax-1, 1:jmax) )
xy_SinLat = sin( xy_Lat )
xy_CosLat = cos( xy_Lat )
! コリオリパラメータの計算
! Calculate Coriolis parameter
!
allocate( xy_Cori (0:imax-1, 1:jmax) )
xy_Cori = 2.0_DP * Omega * xy_SinLat
! 静水圧の式の係数 $ \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_TInpCoefA (1:kmax) )
allocate( z_TInpCoefB (1:kmax) )
allocate( z_TInpCoefK (1:kmax) )
do k = 1, kmax
z_TInpCoefK(k) = &
& ( r_Sigma(k-1) * z_HydroAlpha(k) &
& + r_Sigma(k ) * z_HydroBeta (k) ) &
& / z_DelSigma(k)
enddo
z_TInpCoefA = 0.
do k = 2, kmax
z_TInpCoefA(k) = &
& z_HydroAlpha(k) &
& * ( 1.0_DP - (z_Sigma(k) / z_Sigma(k-1)) ** Kappa ) ** ( -1 )
end do
z_TInpCoefB = 0.
do k = 1, kmax - 1
z_TInpCoefB(k) = &
& z_HydroBeta(k) &
& * ( (z_Sigma(k) / z_Sigma(k+1) ) ** Kappa - 1.0_DP ) ** ( -1 )
end do
! 基準温度 (半整数レベル) の計算
! Calculate reference temperature on half levels
!
allocate( z_RefTemp(1:kmax) )
allocate( r_RefTemp(0:kmax) )
z_RefTemp = RefTemp
r_RefTemp(0) = 0.
r_RefTemp(kmax) = 0.
do k = 1, kmax - 1
r_RefTemp(k) = &
& z_TInpCoefA(k+1) * z_RefTemp(k+1) &
& + z_TInpCoefB( k ) * z_RefTemp( k )
enddo
! ルジャンドル陪関数の固有値の設定
! Set eigen value of Associated Legendre Function
!
allocate( w_LaplaEigVal(lmax) )
w_LaplaEigVal(:) = rn(:,1) / RPlanet**2
! 水平拡散係数の設定
! Configure coefficient of horizontal diffusion
!
! 粘性係数の計算 (最大波数の e-folding time が HDEFoldTime となるように)
! Calculate viscosity coefficient
!
HDEFoldTime = DCCalConvertByUnit( HDEFoldTimeValue, HDEFoldTimeUnit, 'sec' ) ! (in)
if ( HDEFoldTimeValue > 0.0_DP ) then
VisCoef = ( (nmax*(nmax+1)) / RPlanet**2 )**(-HDOrder / 2) &
& / HDEFoldTime
else
VisCoef = 0.0_DP
end if
w_HDifCoefH = - VisCoef * ( ( - w_LaplaEigVal )**( HDOrder / 2 ) )
w_HDifCoefM = w_HDifCoefH - VisCoef * ( - ( 2.0_DP / RPlanet**2 )**( HDOrder / 2 ) )
! スポンジ層の減衰係数の設定
! Set damping coefficient in a sponge layer
!
if ( SLEFoldTimeValue <= 0.0_DP ) then
call MessageNotify( 'E', module_name, &
& 'SLEFoldTimeValue must be greater than zero, SLEFoldTimeValue=<%f>.', &
& d = (/ SLEFoldTimeValue /) )
end if
if ( ( SLNumLayer <= 0 ) .or. ( SLNumLayer > kmax ) ) then
call MessageNotify( 'E', module_name, &
& 'SLNumLayer must be greater than zero and less than/equal to kmax, SLNumLayer=<%d>.', &
& i = (/ SLNumLayer /) )
end if
! スポンジ層の減衰係数の計算
! Calculate damping coefficient for momentum in a sponge layer
!
SLEFoldTime = DCCalConvertByUnit( SLEFoldTimeValue, SLEFoldTimeUnit, 'sec' ) ! (in)
if ( FlagSpongeLayer ) then
! Sponge layer is not applied for model layers, k = 1, kmax-SLNumLayers.
!
do k = 1, kmax-SLNumLayer
wz_SpongeLayerCoefM(:,k) = 0.0_DP
end do
do k = kmax-SLNumLayer+1, kmax
wz_SpongeLayerCoefM(:,k) = &
& 1.0_DP / ( SLEFoldTime * ( z_Sigma(k) / z_Sigma(kmax) )**SLOrder )
end do
if ( .not. FlagSpongeLayerforZonalMean ) then
! Sponge layer is not applied for zonal mean component.
! i.e., sponge layer is applied for zonal wave component only.
!
do k = kmax-SLNumLayer+1, kmax
do l = 1, lmax
a_DegOrd = nm_l( l )
if ( a_DegOrd(2) == 0 ) then
wz_SpongeLayerCoefM(l,k) = 0.0_DP
end if
end do
end do
end if
if ( FlagSpongeLayerforHeat ) then
wz_SpongeLayerCoefH = wz_SpongeLayerCoefM
! Sponge layer is not applied for zonal mean component.
! i.e., sponge layer is applied for zonal wave component only.
!
do k = kmax-SLNumLayer+1, kmax
do l = 1, lmax
a_DegOrd = nm_l( l )
if ( a_DegOrd(2) == 0 ) then
wz_SpongeLayerCoefH(l,k) = 0.0_DP
end if
end do
end do
else
wz_SpongeLayerCoefH = 0.0_DP
end if
else
! Sponge layer is not applied.
!
wz_SpongeLayerCoefM = 0.0_DP
wz_SpongeLayerCoefH = 0.0_DP
end if
! 運動量の水平拡散係数とスポンジ層減衰係数の和
! Damping coefficients by horizonatl diffusion and a sponge layer
!
allocate( wz_DisCoefM(lmax, 1:kmax) )
allocate( wz_DisCoefH(lmax, 1:kmax) )
allocate( w_DisCoefQ (lmax) )
do k = 1, kmax
wz_DisCoefM(:,k) = w_HDifCoefM - wz_SpongeLayerCoefM(:,k)
wz_DisCoefH(:,k) = w_HDifCoefH - wz_SpongeLayerCoefH(:,k)
end do
w_DisCoefQ = w_HDifCoefH
! Calculation of divergence damping period
!
DivDampPeriod = DCCalConvertByUnit( DivDampPeriodValue, DivDampPeriodUnit, 'sec' )
if ( DivDampTime0 > TimeN ) then
call MessageNotify( 'E', module_name, &
& 'DivDampTime0, %f, must be less than or equal to TimeN, %f.', &
& d = (/ DivDampTime0, TimeN /) )
end if
! ヒストリデータ出力のためのへの変数登録
! Register of variables for history data output
!
call HistoryAutoAddVariable( 'DUDtDyn', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'dynamical tendency of zonal wind', 'm s-2' )
call HistoryAutoAddVariable( 'DVDtDyn', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'dynamical tendency of meridional wind', 'm s-2' )
call HistoryAutoAddVariable( 'DTempDtDyn', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'dynamical tendency of temperature', 'K s-1' )
do n = 1, ncmax
call HistoryAutoAddVariable( 'D'//trim(a_QMixName(n))//'DtDyn', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'dynamical tendency of water vapor', 's-1' )
end do
call HistoryAutoAddVariable( 'DPsDtDyn', &
& (/ AxNameX, AxNameY, AxNameT /), &
& 'dynamical tendency of surface pressure', 'Ps s-1' )
call HistoryAutoAddVariable( 'OMG', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'vertical velocity in pressure coordinate (omega, DP/Dt)', 'Pa s-1' )
call HistoryAutoAddVariable( 'SigDot', &
& (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
& 'sigma-vertical velocity', '1 s-1' )
call HistoryAutoAddVariable( 'DPiDt', &
& (/ AxNameX, AxNameY, AxNameT /), &
& 'Pi (log Ps) tendency)', 'Pa s-1' )
call HistoryAutoAddVariable( 'Vor', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'vorticity', 's-1' )
call HistoryAutoAddVariable( 'Div', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'divergence', 's-1' )
call HistoryAutoAddVariable( 'Mass', &
& (/ AxNameX, AxNameY, AxNameT /), &
& 'mass', 'kg' )
call HistoryAutoAddVariable( 'KinEngy', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'kinetic energy', 'J m-2' )
call HistoryAutoAddVariable( 'IntEngy', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'internal energy', 'J m-2' )
call HistoryAutoAddVariable( 'PotEngy', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'potential energy', 'J m-2' )
call HistoryAutoAddVariable( 'LatEngy', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'latent energy', 'J m-2' )
call HistoryAutoAddVariable( 'TotEngy', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'total energy', 'J m-2' )
call HistoryAutoAddVariable( 'Enstro', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'enstrophy', 'kg m-2 s-2' )
! 空間平均された診断量
! Spatially averaged diagnostic variables
!
call HistoryAutoAddVariable( 'TotalMass', &
& (/ AxNameT /), &
& 'total mass', 'kg' )
call HistoryAutoAddVariable( 'TotalKinEngy', &
& (/AxNameT /), &
& 'total kinetic energy', 'J' )
call HistoryAutoAddVariable( 'TotalIntEngy', &
& (/ AxNameT /), &
& 'total internal energy', 'J' )
call HistoryAutoAddVariable( 'TotalPotEngy', &
& (/ AxNameT /), &
& 'total potential energy', 'J' )
call HistoryAutoAddVariable( 'TotalLatEngy', &
& (/ AxNameT /), &
& 'total latent energy', 'J' )
call HistoryAutoAddVariable( 'TotalTotEngy', &
& (/ AxNameT /), &
& 'total total energy', 'J' )
call HistoryAutoAddVariable( 'TotalEnstro', &
& (/ AxNameT /), &
& 'total enstrophy', 'kg m-2 s-2' )
call HistoryAutoAddVariable( 'MeanMass', &
& (/ AxNameT /), &
& 'mean mass', 'kg' )
call HistoryAutoAddVariable( 'MeanKinEngy', &
& (/AxNameT /), &
& 'mean kinetic energy', 'J m-2' )
call HistoryAutoAddVariable( 'MeanIntEngy', &
& (/ AxNameT /), &
& 'mean internal energy', 'J m-2' )
call HistoryAutoAddVariable( 'MeanPotEngy', &
& (/ AxNameT /), &
& 'mean potential energy', 'J m-2' )
call HistoryAutoAddVariable( 'MeanLatEngy', &
& (/ AxNameT /), &
& 'mean latent energy', 'J m-2' )
call HistoryAutoAddVariable( 'MeanTotEngy', &
& (/ AxNameT /), &
& 'mean total energy', 'J m-2' )
call HistoryAutoAddVariable( 'MeanEnstro', &
& (/ AxNameT /), &
& 'mean enstrophy', 'kg m-2 s-2' )
!!$ ! テスト計算出力?
!!$ ! Output for test variable?
!!$ !
!!$ call HistoryAutoAddVariable( 'AdvTestStreamFunc', &
!!$ & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
!!$ & 'mass stream function', 'kg s-1' )
!
! 水平スペクトル解析
! horizontal spectral analysis
!
call DynamicsDiagSpAnalysisInit
! Initialization of modules used in this module
!
! 質量の補正
! Mass fixer
!
call MassFixerInit
if ( FlagAdvTest ) then
!
! Utility module for advection test
!
call AdvTestInit
end if
if ( FlagSLTT ) then
!!$ select case ( IDTTMethod )
!!$ case ( IDTTMethodSp )
!!$ case ( IDTTMethodSL )
! セミラグランジュ法による物質移流
! Semi-Lagrangian method for tracer transport
!
call SLTTInit
!!$ case ( IDTTMethodFV )
!!$ ! 物質移流 (有限体積法)
!!$ ! Tracer Transport (finite volume method)
!!$ !
!!$ call TTFVInit
!!$ end select
end if
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' TimeIntegScheme = %c', c1 = trim( TimeIntegScheme ) )
call MessageNotify( 'M', module_name, ' HDEFoldTime = %f [%c]', &
& d = (/ HDEFoldTimeValue /), c1 = trim(HDEFoldTimeUnit) )
call MessageNotify( 'M', module_name, ' HDOrder = %d', i = (/ HDOrder /) )
call MessageNotify( 'M', module_name, ' VisCoef = %f', d = (/ VisCoef /) )
call MessageNotify( 'M', module_name, ' FlagSpongeLayer = %b', l = (/ FlagSpongeLayer /) )
call MessageNotify( 'M', module_name, ' FlagSpongeLayerforZonalMean = %b', l = (/ FlagSpongeLayerforZonalMean /) )
call MessageNotify( 'M', module_name, ' FlagSpongeLayerforHeat = %b', l = (/ FlagSpongeLayerforHeat /) )
call MessageNotify( 'M', module_name, ' SLEFoldTime = %f [%c]', &
& d = (/ SLEFoldTimeValue /), c1 = trim(SLEFoldTimeUnit) )
call MessageNotify( 'M', module_name, ' SLOrder = %d', i = (/ SLOrder /) )
call MessageNotify( 'M', module_name, ' SLNumLayer = %d', i = (/ SLNumLayer /) )
call MessageNotify( 'M', module_name, ' RefTemp = %f', d = (/ RefTemp /) )
call MessageNotify( 'M', module_name, ' FlagDivDamp = %b', l = (/ FlagDivDamp /) )
call MessageNotify( 'M', module_name, ' DivDampPeriod = %f', d = (/ DivDampPeriod /) )
call MessageNotify( 'M', module_name, ' DivDampTime0 = %f', d = (/ DivDampTime0 /) )
call MessageNotify( 'M', module_name, ' FlagTotMassFix = %b', l = (/ FlagTotMassFix /) )
call MessageNotify( 'M', module_name, ' FlagMassHorDifCor = %b', l = (/ FlagMassHorDifCor /) )
call MessageNotify( 'M', module_name, ' FlagSLTT = %b', l = (/ FlagSLTT /) )
!!$ call MessageNotify( 'M', module_name, ' TracerTransportMethod = %c', c1 = trim( TracerTransportMethod ) )
call MessageNotify( 'M', module_name, ' FlagAdvTest = %b', l = (/ FlagAdvTest /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
dynamics_hspl_vas83_inited = .true.
end subroutine DynamicsHSplVAS83Init
!--------------------------------------------------------------------------------------
subroutine DynamicsHSplVAS83Finalize
!
! モジュール内部の変数の割り付け解除を行います.
!
! Deallocate variables in this module.
!
! 宣言文 ; Declaration statements
!
implicit none
! 実行文 ; Executable statement
!
if ( .not. dynamics_hspl_vas83_inited ) return
! デフォルト値へ戻す
! Return to default values
!
DelTimeSave = - 1.0_DP
! 割り付け解除
! Deallocation
!
if ( allocated( xy_SinLat ) ) deallocate( xy_SinLat )
if ( allocated( xy_CosLat ) ) deallocate( xy_CosLat )
if ( allocated( xy_Cori ) ) deallocate( xy_Cori )
if ( allocated( z_HydroAlpha ) ) deallocate( z_HydroAlpha )
if ( allocated( z_HydroBeta ) ) deallocate( z_HydroBeta )
if ( allocated( z_TInpCoefA ) ) deallocate( z_TInpCoefA )
if ( allocated( z_TInpCoefB ) ) deallocate( z_TInpCoefB )
if ( allocated( z_TInpCoefK ) ) deallocate( z_TInpCoefK )
if ( allocated( z_RefTemp ) ) deallocate( z_RefTemp )
if ( allocated( r_RefTemp ) ) deallocate( r_RefTemp )
if ( allocated( w_LaplaEigVal ) ) deallocate( w_LaplaEigVal )
if ( allocated( wz_DisCoefM ) ) deallocate( wz_DisCoefM )
if ( allocated( wz_DisCoefH ) ) deallocate( wz_DisCoefH )
if ( allocated( w_DisCoefQ ) ) deallocate( w_DisCoefQ )
if ( allocated( z_siMtxC ) ) deallocate( z_siMtxC )
if ( allocated( z_siMtxG ) ) deallocate( z_siMtxG )
if ( allocated( zz_siMtxH ) ) deallocate( zz_siMtxH )
if ( allocated( zz_siMtxDiH ) ) deallocate( zz_siMtxDiH )
if ( allocated( wzz_siMtxWDiH ) ) deallocate( wzz_siMtxWDiH )
if ( allocated( zz_siMtxGCt ) ) deallocate( zz_siMtxGCt )
if ( allocated( zz_siMtxW ) ) deallocate( zz_siMtxW )
if ( allocated( zz_siMtxQ ) ) deallocate( zz_siMtxQ )
if ( allocated( zz_siMtxS ) ) deallocate( zz_siMtxS )
if ( allocated( zz_siMtxR ) ) deallocate( zz_siMtxR )
!!$ if ( allocated(nmo ) ) deallocate( nmo )
!!$ if ( allocated(wzz_siMtxM ) ) deallocate( wzz_siMtxM )
!!$ if ( allocated(z_siMtxPivWork) ) deallocate( z_siMtxPivWork )
if ( allocated(wzz_siMtxLU ) ) deallocate( wzz_siMtxLU )
if ( allocated(wz_siMtxPiv ) ) deallocate( wz_siMtxPiv )
dynamics_hspl_vas83_inited = .false.
end subroutine DynamicsHSplVAS83Finalize
!-------------------------------------------------------------------
! Memo: YOT (2009/10/04)
! The routine, SemiImplMatrix_old, is old version of SemiImplMatrix.
! This is not deleted now, since this may be valuable as a reference in
! optimizing the code.
!
!!$
!!$ subroutine SemiImplMatrix_old
!!$ !
!!$ ! TimeIntegration で使用する係数の設定を行います.
!!$ ! 初回および $ \Delta t $ が変更された場合以外は,
!!$ ! 前回に設定した値をそのまま用います.
!!$ !
!!$ ! Configure coefficients for "TimeIntegration".
!!$ ! Setting values that are set last time are used,
!!$ ! except when first time and Δt are changed.
!!$ !
!!$
!!$ ! モジュール引用 ; USE statements
!!$ !
!!$
!!$ ! 物理定数設定
!!$ ! Physical constants settings
!!$ !
!!$ use constants, only: &
!!$ & RPlanet, &
!!$ ! $ a $ [m].
!!$ ! 惑星半径.
!!$ ! Radius of planet
!!$ & CpDry
!!$ ! $ C_p $ [J kg-1 K-1].
!!$ ! 乾燥大気の定圧比熱.
!!$ ! Specific heat of air at constant pressure
!!$
!!$ ! 座標データ設定
!!$ ! Axes data settings
!!$ !
!!$ use axesset, only: &
!!$ & r_Sigma, & ! $ \sigma $ レベル (半整数).
!!$ ! Half $ \sigma $ level
!!$ & 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 次方程式を解くための関数 (spml 同梱モジュール)
!!$ ! Functions to solve linear simultaneous equation by LU decomposition
!!$ ! (a module included in spml)
!!$ !
!!$ 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.
!!$ !
!!$ 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 ( DelTimeSave == DelTime ) return
!!$ DelTimeSave = DelTime
!!$
!!$ call DbgMessage( '%c: %c: (DelTime=%f [sec]) coefficients for "TimeIntegration" is generated. ', &
!!$ & c1 = module_name, c2 = 'SemiImplMatrix', 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_TInpCoefK * z_RefTemp
!!$
!!$ 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.
!!$ 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.
!!$ 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.
!!$ do k = 1, kmax
!!$ zz_siMtxQ(k,k) = ( r_RefTemp(k-1) - z_RefTemp(k) ) / z_DelSigma(k)
!!$ enddo
!!$ do k = 1, kmax-1
!!$ zz_siMtxQ(k,k+1) = ( z_RefTemp(k) - r_RefTemp(k) ) / z_DelSigma(k)
!!$ enddo
!!$
!!$ zz_siMtxQS = 0.
!!$ zz_siMtxQS = matmul(zz_siMtxQ, zz_siMtxS)
!!$
!!$ zz_siMtxR = 0.
!!$ do k = 1, kmax
!!$ do l = k, kmax
!!$ zz_siMtxR(k,l) = &
!!$ & - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(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_RefTemp(k)
!!$ enddo
!!$ enddo
!!$
!!$ zz_siMtxH = 0.
!!$ zz_siMtxH = zz_siMtxQS - zz_siMtxR
!!$
!!$ zz_siMtxWH = 0.
!!$ zz_siMtxWH = matmul(zz_siMtxW, zz_siMtxH)
!!$
!!$ if ( .not. allocated(nmo ) ) allocate( nmo (1:2, 0:nmax, 0:nmax) )
!!$ if ( .not. allocated(wzz_siMtxM ) ) allocate( wzz_siMtxM ((nmax+1)**2, 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
!!$
!!$ ! スペクトル添字順序の取り出し
!!$ ! 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
!!$ wzz_siMtxM ( nm,k,kk ) = &
!!$ & - DelTime**2 * wz_LaplaEigVal(nm,1) &
!!$ & * ( zz_siMtxWH( k,kk ) &
!!$ & + zz_siMtxGCt( k,kk ) &
!!$ & * ( 1. - 2. * DelTime * wz_HDifCoefH(nm,1) ) )
!!$ if ( k == kk ) then
!!$ wzz_siMtxM ( nm,k,kk ) = &
!!$ & wzz_siMtxM ( nm,k,kk ) &
!!$ & + ( 1. - 2. * DelTime * wz_HDifCoefM(nm,1) ) &
!!$ & * ( 1. - 2. * DelTime * wz_HDifCoefH(nm,1) )
!!$ endif
!!$ end do
!!$ end do
!!$
!!$ ! LU 行列計算
!!$ ! LU matrix calculation
!!$ !
!!$ call LUDecomp( &
!!$ & wzz_siMtxM(nm,:,:), & ! (inout)
!!$ & z_siMtxPivWork ) ! (out)
!!$
!!$ ! ダミー値の代入. (位置 kmax はまだ未定義なため).
!!$ ! Dummy value is subtituted. (Because position kmax is undefined yet).
!!$ !
!!$ z_siMtxPivWork(kmax) = 0
!!$
!!$ ! 配列の詰め替え
!!$ ! 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 ) = wzz_siMtxM ( nm,k,kk )
!!$ wzz_siMtxLU ( nmo(2,m,l),k,kk ) = wzz_siMtxM ( nm,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 ) = wzz_siMtxM ( nm,k,kk )
!!$ end do
!!$ wz_siMtxPiv ( nm,k ) = z_siMtxPivWork ( k )
!!$ end do
!!$ end do
!!$
!!$ end subroutine SemiImplMatrix_old
!!$
!--------------------------------------------------------------------------------------
subroutine SemiImplMatrix
!
! TimeIntegration で使用する係数の設定を行います.
! 初回および $ \Delta t $ が変更された場合以外は,
! 前回に設定した値をそのまま用います.
!
! Configure coefficients for "TimeIntegration".
! Setting values that are set last time are used,
! except when first time and $ \Delta t $ are changed.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: &
& RPlanet, &
! $ a $ [m].
! 惑星半径.
! Radius of planet
& CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax ! 鉛直層数.
! Number of vertical level
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& r_Sigma, & ! $ \sigma $ レベル (半整数).
! Half $ \sigma $ level
& 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)
!
#ifdef AXISYMMETRY
use wa_zonal_module, only: l_nm
#elif LIB_MPI
use ua_mpi_module, only: l_nm => lu_nm
#else
use wa_module, only: l_nm
#endif
! LU 分解法により連立 1 次方程式を解くための関数 (spml 同梱モジュール)
! Functions to solve linear simultaneous equation by LU decomposition
! (a module included in spml)
!
use lumatrix, only: LUDecomp
! デバッグ用ユーティリティ
! Utilities for debug
!
use dc_trace, only: DbgMessage, BeginSub, EndSub, Debug
! 宣言文 ; Declaration statements
!
implicit none
! TimeIntegration 等で使用する係数の設定のための作業変数
! Work variable for coefficients for "TimeIntegration", etc.
!
integer:: k, l, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n
! 実行文 ; Executable statement
!
! $ \Delta t $ [s] のチェック・保存
! Check and save $ \Delta t $ [s]
!
if ( DelTimeSave == DelTime ) return
DelTimeSave = DelTime
call DbgMessage( '%c: %c: (DelTime=%f [sec]) coefficients for "TimeIntegration" is generated. ', &
& c1 = module_name, c2 = 'SemiImplMatrix', d = (/ DelTime /) )
! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
if ( .not. allocated( z_siMtxC ) ) allocate( z_siMtxC (1:kmax) )
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( wz_siMtxDi ) ) allocate( wz_siMtxDi (lmax,1:kmax) )
if ( .not. allocated( zz_siMtxDiH ) ) allocate( zz_siMtxDiH (1:kmax, 1:kmax) )
if ( .not. allocated( wzz_siMtxWDiH ) ) allocate( wzz_siMtxWDiH(lmax,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_siMtxR ) ) allocate( zz_siMtxR (1:kmax, 1:kmax) )
z_siMtxC = z_DelSigma
z_siMtxG = CpDry * z_TInpCoefK * z_RefTemp
do k = 1, kmax
do l = 1, kmax
zz_siMtxGCt(k,l) = z_siMtxG(k) * z_siMtxC(l)
end do
end do
zz_siMtxW = 0.
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.
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.
do k = 1, kmax
zz_siMtxQ(k,k) = ( r_RefTemp(k-1) - z_RefTemp(k) ) / z_DelSigma(k)
enddo
do k = 1, kmax-1
zz_siMtxQ(k,k+1) = ( z_RefTemp(k) - r_RefTemp(k) ) / z_DelSigma(k)
enddo
zz_siMtxR = 0.
do k = 1, kmax
do l = k, kmax
zz_siMtxR(k,l) = &
& - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(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_RefTemp(k)
enddo
enddo
zz_siMtxH = matmul(zz_siMtxQ, zz_siMtxS) - zz_siMtxR
! Check of coefficients for horizontal diffusion and sponge layer
! A threshold value used below is arbitrary.
!
do n = 1, lmax
do k = 1, kmax
if ( abs( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM(n,k) ) < 1.0e-10_DP ) then
call MessageNotify( 'E', module_name, &
& 'Dissipation coefficient for momentum is inappropriate.' )
end if
if ( abs( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH(n,k) ) < 1.0e-10_DP ) then
call MessageNotify( 'E', module_name, &
& 'Dissipation coefficient for heat is inappropriate.' )
end if
end do
if ( abs( 1.0_DP - 2.0_DP * DelTime * w_DisCoefQ(n) ) < 1.0e-10_DP ) then
call MessageNotify( 'E', module_name, &
& 'Dissipation coefficient for composition is inappropriate.' )
end if
end do
wz_siMtxDi = 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH )
do n = 1, lmax
! NOTE: wz_siMtiDi is a diagonal matrix.
!
do k = 1, kmax
do kk = 1, kmax
zz_siMtxDiH(k,kk) = wz_siMtxDi(n,k) * zz_siMtxH(k,kk)
end do
end do
zz_siMtxDiH = matmul( zz_siMtxW, zz_siMtxDiH )
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxWDiH(n,k,kk) = zz_siMtxDiH(k,kk)
end do
end do
end do
if ( .not. allocated(wzz_siMtxLU) ) allocate( wzz_siMtxLU(lmax, 1:kmax, 1:kmax) )
if ( .not. allocated(wz_siMtxPiv) ) allocate( wz_siMtxPiv(lmax, 1:kmax) )
! 行列 $ \underline{M} $ の計算
! Calculate matrix $ \underline{M} $
!
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxLU ( :,k,kk ) = &
& - DelTime**2 * w_LaplaEigVal(:) &
& * ( wzz_siMtxWDiH(:,k,kk) + zz_siMtxGCt(k,kk) )
if ( k == kk ) then
wzz_siMtxLU ( :,k,kk ) = &
& wzz_siMtxLU ( :,k,kk ) &
& + ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM (:,k) )
endif
end do
end do
! LU 行列計算
! LU matrix calculation
!
call LUDecomp( &
& wzz_siMtxLU, & ! (inout)
& wz_siMtxPiv ) ! (out)
end subroutine SemiImplMatrix
!--------------------------------------------------------------------------------------
end module dynamics_hspl_vas83