3. 力学過程

3.1 力学過程の概要

本章では力学過程のコードの構造と計算手順を説明する. ここでいう力学過程と はスペクトルで計算される部分のことを指し, 移流項や水平拡散項, 圧力勾配項 などである. あるいは物理過程以外と言っても良いかもしれない. スペクトル とは球面調和関数展開の係数のことであり, スペクトルへの変換とは球関数展開 に他ならない.

実際の時間積分は力学過程の中で行われる. 物理過程で計算された時間変化 項は力学過程に引き渡され, 力学過程の時間変化項に加えられ, セミインプ リシット法による時間積分に組み込まれる. ただし, 降水過程などの調節過程は時 間積分後に適用される.

3.2 力学過程 Dynamics の流れ

Fig.3.1〜 Fig.3.3 に力学過程サブルーチン Dynamics の流れを示す. 四角枠内の太文字はサブルーチン名を表している. Dynamics とサブルーチン間にやりとりされる変数が示されている. 四角枠内が空白であるものは, 簡略に記述できるため Dynamics サブルーチン本体に直接コードが記述されているものである. これらのサブルーチンは, モジュール dynamics_hspl_vas83 内に ひとまとめにされている. DynamicsInitおよびSemiImplMatrix 内で設定されたパラメータや配列はモジュール内の共有変数として保存され, 各サブルーチンはその変数を参照する仕組みとなっているため, 原則的に Dynamics と各サブルーチンとでやり取りされるのは 時間変化する変数のみとなる.

Dynamics で呼ばれている NonLinearOnGrid は時間変化率の非重力波(NG)項を計算する. TimeIntegration はセミインプリシット時間積分 (3.3.1節参照)を行う. 以上の2つが Dynamics の主な下請ルーチンである.

NonLinearOnGrid 以前に呼ばれるサブルーチンは 必要な定数等の準備にあてられる. SemiImplMatrix までは原則的に時間に依存しない. TimeIntegration の前後で格子点値とスペクトル値の変換を行う. その後, 水平粘性に失われた運動エネルギーを補償する 摩擦熱を行い, Dynamics 外に出ない量の出力(DiagOutput)を行う. DiagOutput では診断量の計算を行うが, これは モニターのためであり, 時間積分には影響しない.

図 3.1: 力学過程 Dynamics の流れ1. 四角枠内の太文字はサブルーチン名であり, Dynamics が呼ぶサブルーチンと引数を示している. 引数は, 特に注記されていないものは格子点値である. 四角枠内が空白であるものは, 簡略に記述できるため Dynamics サブルーチン本体に直接コードが記述されているものである. ここでは Dynamics が引数として値を受け取り, 初期設定 (DynamicsInit, SemiImplMatrix) および 非線形項の計算 (NonLinearOnGrid) の結果得られた 値からスペクトル値で表現された時間変化を求める ところまでが記されている. これ以後の処理については Fig.3.2を参照されたい.

図 3.2: 力学過程 Dynamics の流れ2. Fig.3.1の続きである. 図の見方についてはFig.3.1を参照せよ. ここでは格子点値をスペクトル値に変換し, それらを用いて, セミインプリシット時間積分 TimeIntegration を行っている. これ以後の処理については Fig.3.3を参照されたい.

図 3.3: 力学過程 Dynamics の流れ3. Fig.3.2の続きである. 図の見方についてはFig.3.1を参照せよ. ここではスペクトル値を格子点値に戻し, 摩擦熱補正をかけた後, ヒストリデータを出力し, 最後に $t+\Delta t$の値をサブルーチンの外に返している.

3.3 セミインプリシット時間積分

注: 本節は書き換え中であるため, 支配方程式系とその離散化の文章とは対応していない.

3.3.1 セミインプリシット時間積分の概要

セミインプリシット法については, 支配方程式系とその離散化ドキュメント の 3.5 節に解説があるので 参照のこと. 時刻 $t-\Delta t$ から $t+\Delta t$ へのセミインプリシットの 計算手順は以下のようにまとめられる.

1. 発散項 $\Dvect{D}$ のみに関する単一のセミインプリシット方程式,

   $$\tilde{\underline{M}}^{m}_{n}\, \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} = \tilde{\Dvect{f}}^{m}_{n}$$ (3.1)

   
   
   を解く. ここで
   $$\displaystyle \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} \equiv \f... ...e{\Dvect{D}}^{m,t+\Delta t}_{n} + \tilde{\Dvect{D}}^{m,t-\Delta t}_{n} \right)$$
   である.

2. $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ と時刻 $t$ における時間変化の非重力波 (NG)項から各変数の $t+\Delta t$ での値を求める.

以上の手順の説明は 3.3.2 節で行う.

3.3.2 セミインプリシット時間積分の詳細

3.3.1節で述べたように, セミインプリシット時間積分 は2つのステップに分けられる. これを詳しく書くと以下のようになる.

● 第1段階: $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ の計算

(3.1) を解くのが第1段階 である. ここでは支配方程式系とその離散化ドキュメント の 3.5 節と同様, 太字で鉛直に離散化したベクトルを表 し, 下線で行列を表す. $\Dvect{D}$ は発散で,

$$\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} \equiv \frac{1}{2} (\tilde{\Dvect{D}}^{m,t+\Delta t}_{n} + \tilde{\Dvect{D}}^{m,t-\Delta t}_{n})$$ (3.2)

である. 行列 $\tilde{\underline{M}}^{m}_{n}$ は

$$\begin{align*}\begin{split}\tilde{\underline{M}}^{m}_{n} &\equiv ( 1 -2\Delta t ... ...vect{C}^{T} \right\} \left( - \frac{n(n+1)}{a^2} \right) \end{split}\end{align*}$$ (3.3)

$$\begin{align*}\begin{split}\tilde{\Dvect{f}}^{m}_{n} & = ( 1 -2\Delta t \tilde{\... ...tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} \right] \right\} \end{split}\end{align*}$$ (3.4)

である.

$\tilde{\underline{M}}^{m}_{n}$ はサブルーチン SemiImplMatrix で計算され, LU分解される. $\tilde{\underline{M}}^{m}_{n}$ は時間刻みの大きさが変らない限り設定し直す必要がない. $\tilde{\Dvect{f}}^{m}_{n}$ の計算と (3.1) を解く作業 はサブルーチン TimeIntegration で行われる.

● 第2段階: 時間積分
時間積分も TimeIntegration が行う. 第1段階で $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ が 求まったので発散の時間積分は容易である:

$$\tilde{\Dvect{D}}^{m,t+\Delta t}_{n} = 2 \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} - \tilde{\Dvect{D}}^{m,t-\Delta t}_{n} .$$ (3.5)

他の物理量についても, $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ と時刻 $t$ における時間変化率の非重力波(NG)項から, $t+\Delta t$ での値が求まる:

$$\tilde{\cal X}^{m,t+\Delta t}_n = (\tilde{\gamma}_{{\cal X},n}^{m... ...{{\cal X},n}^{m} \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} \right] \right\} \, .$$ (3.6)

ここで,

$$\begin{array}{cccccc} \displaystyle \tilde{\cal X}^{m}_{n} & ... ...tyle -\underline{h}, & \displaystyle \underline{0}. \end{array}$$ (3.7)

以下の各節で, 具体的な計算手順とプログラムソースとの対応を述べる.

3.4 モデル変数表

ここでは, これ以後の節で述べられるモデル変数の表記法を述べる. モデル変数は    $$ のようなフォントで示され, その入出力に関連して表 3.1 のように分類される.

表 3.1: 数式と対応させて書く場合の モデル変数の表記法. 本文中で定数の下に付した文字列は, 値を設定したサブルーチン名を示す. 緑で示される定数は, 物理定数管理, 座標データ管理モジュールなどから参照するか, もしくはそれらの値を用いて DynamicsInit や SemiImplMatrix で時間ループの初回に設定される.

解説	例
height 1pt 当該サブルーチンからの出力変数 (入力に上書きする場合を含む) の出力時	$^{\mbox{$\bullet$}}$
当該関数からの返り値	〃
当該サブルーチン内の作業変数または, 入(出)力変数が入(出)力時の値と違う 場合	$$
当該サブルーチンへの入力変数	$$
当該関数への入力変数	〃
当該サブルーチンへの入力に上書きして出力される変数の入力時	$^{\mbox{$\bullet$}}$
定数	$$
関数	$$

3.5 非重力波項の計算 (サブルーチン NonLinearOnGrid)

3.5.1 NonLinearOnGrid前の準備

NonLinearOnGridの直前で計算される, 地表面気圧の空間変化と渦度発散について記す. 灰色文字で示されるのは SPMODEL ライブラリ (Takehiro et al., 2006) から提供される 関数である. この節では変数の時間ステップは全て $t$ であるため, 簡単化のため $(t)$ の表記を省略する. 経度, 緯度, 鉛直方向添え字 $i,j,k$ の表記も省略する.

$$\underbrace{\left(\nabla \pi \right)_{\lambda}}_{ \mbox{{\cmssbx xy\_GradLambdaPiN}}^{\mbox{$\bullet$}} } = \underbrace{\Dinv{a}}_{ \mbox{{\cmssbx\textcolor{PineGreen}{RPl... ...lor{Gray}{\qquad\qquad\qquad\qquad\qquad xy\_GradLambda\_w( w\_xy( ... ) )}}} }$$ (3.8)

$$\underbrace{\left(\nabla \pi \right)_{\mu}}_{ \mbox{{\cmssbx xy\_GradMuPiN}}^{\mbox{$\bullet$}} } = \underbrace{\Dinv{a}}_{ \mbox{{\cmssbx\textcolor{PineGreen}{RPl... ...x{{\cmssbx\textcolor{Gray}{\qquad\qquad\quad xy\_GradMu\_w( w\_xy( ... ) )}}} }$$ (3.9)

$$\underbrace{\zeta }_{ \mbox{{\cmssbx xyz\_VorN}}^{\mbox{$\bullet$}} } = \underbrace{\Dinv{a}}_{ \mbox{{\cmssbx\textcolor{PineGreen}{RPl... ...!\!\!\!\!\!\!\!} \biggr) }_{ \mbox{{\cmssbx\textcolor{Gray}{xya\_wa( ... )}}} }$$ (3.10)

$$\underbrace{D }_{ \mbox{{\cmssbx xyz\_VorN}}^{\mbox{$\bullet$}} } = \underbrace{\Dinv{a}}_{ \mbox{{\cmssbx\textcolor{PineGreen}{RPl... ...!\!\!\!\!\!\!\!} \biggr) }_{ \mbox{{\cmssbx\textcolor{Gray}{xya\_wa( ... )}}} }$$ (3.11)

3.5.2 NonLinearOnGridでの計算

サブルーチン NonLinearOnGrid では格子点上で $U_A, V_A, H, R, KE,$ $UT', VT',$ $Uq, Vq,$ $\dot{\sigma}, Z$ について計算を行う. この節では時間 $(t)$ と 経度, 緯度方向添え字 $i,j$ の表記を省略する.

地表面気圧の変化に関連する項の計算

$$\underbrace{ Z }_{ \mbox{{\cmssbx xy\_DPiDtN}}^{\mbox{$\bullet$}} } = - \underbrace{ \sum_{k=1}^{K} \overbrace{\Dvect{v}_{k} \cdot \n... ...tcolor{PineGreen}{z\_DelSigma}}} } }_{ \mbox{{\cmssbx xyz\_PiAdvSum$_{k=1}$}} }$$ (3.12)

ここで,

$$\underbrace{\Dvect{v}_{k} \cdot \nabla \pi}_{ \mbox{{\cmssbx xyz\_PiAdv}} } = \overbrace{U_k}^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!... ...1-\mu^2)}{a}\DP{\pi}{\mu}.}^{ \mbox{{\cmssbx\textcolor{blue}{xy\_GradMuPiN}}} }$$ (3.13)

渦度, 発散の変化に関連する項の計算

$$\begin{align*}\begin{split}\overbrace{{U_A}_{,1}}^{\!\!\!\!\!\!\!\! \mbox{{\cmss... ... \DP{\pi}{\lambda} + {\cal F}_{\lambda, K} \cos \varphi, \end{split}\end{align*}$$ (3.14)

$$\begin{align*}\begin{split}\overbrace{{V_A}_{,1}}^{\!\!\!\!\!\!\!\! \mbox{{\cmss... ...}{a} \DP{\pi}{\mu} + {\cal F}_{\varphi, K} \cos \varphi, \end{split}\end{align*}$$ (3.15)

$$\underbrace{ (\mbox{\sl KE})_{k} + \sum_{l=1}^{K}W_{kl} (T_{v,l}-... ...}_{\!\!\!\!\!\! \mbox{\footnotesize {\tt HydroGrid}を呼ぶと $W_{kl}$がかかる} }$$ (3.16)

ここで,

$$\underbrace{U_k}_{ \mbox{{\cmssbx xyz\_UCosLatN}} } = \overbrace{u_k}^{ \mbox{{\cmssbx\textcolor{blue}{xyz\_UN}}} } \... ...}_{\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{xy\_CosLat}}} \!\!\!\!\!\!}$$ (3.17)

$$\underbrace{V_k}_{ \mbox{{\cmssbx xyz\_VCosLatN}} } = \overbrace{v_k}^{ \mbox{{\cmssbx\textcolor{blue}{xyz\_VN}}} } \... ...}_{\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{xy\_CosLat}}} \!\!\!\!\!\!}$$ (3.18)

$$\begin{align*}\begin{split}\underbrace{\dot{\sigma}_{k-1/2}}_{ \mbox{{\cmssbx xy... ...\ \dot{\sigma}_{1/2} &= 0, \\ \dot{\sigma}_{K+1/2} &= 0, \end{split}\end{align*}$$ (3.19)

$$\underbrace{T_{v,k}^{\prime}}_{ \mbox{{\cmssbx xyz\_TempVirEdd}} } = \overbrace{T_{v,k}}^{ \mbox{{\cmss... ...{T}_k}_{\!\!\!\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{z\_RefTemp}}} },$$ (3.20)

$$\underbrace{T_{v,k}}_{ \mbox{{\cmssbx xyz\_TempVir}} } = \underbrace{T_k}_{ \mbox{{\cmssbx\textcolor{blue}{xyz\_TempN}}}... ...\mbox{{\cmssbx\textcolor{blue}{\!\!\!\!\!\!\!\!\!\!\!\!xyz\_QVapN}}} } \bigg\},$$ (3.21)

$$\begin{align*}\begin{split}\underbrace{\hat{\kappa}_k}_{\!\!\!\!\!\!\!\!\!\!\!\!... ...\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!} } { \Delta \sigma_k }, \end{split}\end{align*}$$ (3.22)

$$\begin{align*}\begin{split}\underbrace{\alpha_k}_{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!... ...lor{PineGreen}{Kappa}}} \!\!\!\!\!\!\!\!\!\!\!\!\!}} -1, \end{split}\end{align*}$$ (3.23)

$$\begin{align*}\begin{split}\underbrace{\beta_k}_{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\... ...( \frac{ \sigma_{k+1/2} } { \sigma_k } \biggl)^{\kappa}, \end{split}\end{align*}$$ (3.24)

$$\begin{align*}\begin{split}\underbrace{\kappa}_{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!... ...{PineGreen}{CpDry}}} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}}. \end{split}\end{align*}$$ (3.25)

温度の変化に関連する項の計算

$$^{\mbox{$\bullet$}} \equiv \underbrace{U_k}_{ \mbox{{\cmssbx xyz\_UCosLatN}} } \!\!\!... ...mssbx xyz\_TempEdd}} \ \mbox{\scriptsize\Deqref{code-nongrav:温度の擾乱}参照} }$$ (3.26)

$$^{\mbox{$\bullet$}} \equiv \underbrace{V_k}_{ \mbox{{\cmssbx xyz\_VCosLatN}} } \!\!\!\!\!\!\!\!\!\! T'_k,$$ (3.27)

$$\begin{align*}\begin{split}\underbrace{ H_k }_{ \mbox{{\cmssbx xyz\_TempNonLinea... ...lta \sigma_{K} + T'_{v,K} D_K \Delta \sigma_{K} \right]. \end{split}\end{align*}$$ (3.28)

ここで,

$$\underbrace{T_{k}^{\prime}}_{ \mbox{{\cmssbx xyz\_TempEdd}} } = \overbrace{T_{k}}^{ \mbox{{\cmssbx\textcolor{blue}{xyz\_TempN}}... ...{T}_k}_{\!\!\!\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{z\_RefTemp}}} },$$ (3.29)

$$\underbrace{\hat{T}_{k-1/2}^{\prime}}_{ \mbox{{\cmssbx xyr\_TempE... ...} & \text{($k = 2, \cdots, K$)} \\ 0 , & \text{($k = K+1$)} \end{array} \right.$$ (3.30)

$$\underbrace{\hat{\overline{T}}_{k-1/2}}_{ \begin{array}{r} \mbox{... ..., & \text{($k = 2, \cdots, K$)} \\ 0 , & \text{($k = K+1$)} \end{array} \right.$$ (3.31)

$$\underbrace{a_k}_{\begin{array}{r} \mbox{{\cmssbx\textcolor{PineGreen}{z\_TInpCoefA}}} \\ \mbox{{\footnotesize\tt DynamicsInit}} \end{array}} = \overbrace{\alpha_k}^{ \mbox{{\cmssbx\textcolor{PineGreen}{z\_H... ...eft[ 1- \left( \frac{ \sigma_k }{ \sigma_{k-1} } \right)^{\kappa} \right]^{-1},$$ (3.32)

$$\underbrace{b_k}_{\begin{array}{r} \mbox{{\cmssbx\textcolor{PineGreen}{z\_TInpCoefB}}} \\ \mbox{{\footnotesize\tt DynamicsInit}} \end{array}} = \overbrace{\beta_k}^{ \mbox{{\cmssbx\textcolor{PineGreen}{z\_Hy... ...ft[ \left( \frac{ \sigma_k }{ \sigma_{k+1} } \right)^{\kappa} - 1 \right]^{-1},$$ (3.33)

$$\underbrace{ \dot{\sigma}^{\rm NG}_{k-1/2}}_{ \mbox{{\cmssbx xyr\... ...{l} \cdot \nabla \pi \Delta \sigma_{l}}_{ \mbox{{\cmssbx xyz\_PiAdvSum}} } \, .$$ (3.34)

比湿の変化に関連する項の計算

$$^{\mbox{$\bullet$}} \equiv \underbrace{U_k}_{ \mbox{{\cmssbx xyz\_UCosLatN}} } \!\!\!... ...!\!\!\!\!\! \overbrace{q_{k},}^{ \mbox{{\cmssbx\textcolor{blue}{xyz\_QVapN}}} }$$ (3.35)

$$^{\mbox{$\bullet$}} \equiv \underbrace{V_k}_{ \mbox{{\cmssbx xyz\_VCosLatN}} } \!\!\!\!\!\!\!\! q_k,$$ (3.36)

$$\begin{align*}\begin{split}\overbrace{R_1}^{ \mbox{{\cmssbx xyz\_QVapNonLinearN}... ...Delta \sigma_K} \dot{\sigma}_{K-1/2} ( q_{K-1} - q_K ) . \end{split}\end{align*}$$ (3.37)

3.5.3 スペクトル時間変化項の計算

サブルーチン NonLinearOnGrid から得られた値から $\pi, \zeta, D, T, q$に関するスペクトル時間変化項を計算する. この節でも変数の時間ステップは全て $t$ であるため, 簡単化のため $(t)$ の表記を省略する. なお, 本来 $\widetilde{[ \hspace{0.5cm} ]}^{m}_{n}$と記載されるべき 表記を, 変数名を記載する関係からチルダを省略して $[ \hspace{0.5cm} ]^{m}_{n}$と記載する.

$$\underbrace{ \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} ... ...} } \!\!\!\!\!\!\!\bigr]^{m}_{n}, }^{ \mbox{{\cmssbx\textcolor{Gray}{w\_xy}}} }$$ (3.38)

$$\begin{align*}\begin{split}\overbrace{ \biggl( \DP{\tilde{\zeta}^{m}_{k,n}}{t} \... ...ox{{\cmssbx\textcolor{Gray}{wa\_DivMu\_xya}}} } \biggr), \end{split}\end{align*}$$ (3.39)

$$\begin{align*}\begin{split}\overbrace{ \biggl( \DP{\tilde{D}^{m}_{k,n}}{t} \bigg... ...x\textcolor{Gray}{wa\_Lapla\_wa( wa\_xya( ... ) )}}} } , \end{split}\end{align*}$$ (3.40)

$$\begin{align*}\begin{split}\overbrace{ \biggl( \DP{\tilde{T}^{m}_{k,n}}{t} \bigg... ...m}_{n} }_{ \mbox{{\cmssbx\textcolor{Gray}{wa\_xya}}} } , \end{split}\end{align*}$$ (3.41)

$$\begin{align*}\begin{split}\overbrace{ \biggl( \DP{\tilde{q}^{m}_{k,n}}{t} \bigg... ...m}_{n} }_{ \mbox{{\cmssbx\textcolor{Gray}{wa\_xya}}} } , \end{split}\end{align*}$$ (3.42)

3.6 格子点値からスペクトル値への変換

地形, $t-\Delta t$の予報変数    $,$   $$, および 物理過程の時間変化項    $,$   $$ の格子点値から, スペクトル値を算出する. この節でも3.5.3節と同様に, $\widetilde{[ \hspace{0.5cm} ]}^{m}_{n}$と記載されるべき 表記を, 変数名を記載する関係からチルダを省略して $[ \hspace{0.5cm} ]^{m}_{n}$と記載する.

3.6.1 地形

$$\underbrace{ \Phi^{m}_{s,n} }_{ \mbox{{\cmssbx w\_SurfGeoPot}}^{\... ...!\!\!\!\!\!} \!\! \bigr]^{m}_{n}. }^{ \mbox{{\cmssbx\textcolor{Gray}{w\_xy}}} }$$ (3.43)

3.6.2 $t-\Delta t$の予報変数

この節の変数はすべて時間ステップ $t-\Delta t$ の値である. 簡単化のため $(t-\Delta t)$ の表記を省略する.

$$\begin{align*}\begin{split}\underbrace{ \tilde{\zeta}^{m}_{k,n} }_{ \mbox{{\cmss... ...ox{{\cmssbx\textcolor{Gray}{wa\_DivMu\_xya}}} } \biggr), \end{split}\end{align*}$$ (3.44)

$$\begin{align*}\begin{split}\underbrace{ \tilde{D}^{m}_{k,n} }_{ \mbox{{\cmssbx w... ...x{{\cmssbx\textcolor{Gray}{wa\_DivMu\_xya}}} } \biggr) , \end{split}\end{align*}$$ (3.45)

$$\underbrace{ T^{m}_{k,n} }_{ \mbox{{\cmssbx wz\_TempB}}^{\mbox{$\bullet$}} } = \overbrace{ \bigl[ \!\! \underbrace{ T_{ijk} }_{ \mbox{{\cmssbx... ...TempB}}} } \!\! \bigr]^{m}_{n}, }^{ \mbox{{\cmssbx\textcolor{Gray}{wa\_xya}}} }$$ (3.46)

$$\underbrace{ q^{m}_{k,n} }_{ \mbox{{\cmssbx wz\_QVapB}}^{\mbox{$\bullet$}} } = \overbrace{ \bigl[ \!\! \underbrace{ q_{ijk} }_{ \mbox{{\cmssbx... ...QVapB}}} } \!\! \bigr]^{m}_{n}, }^{ \mbox{{\cmssbx\textcolor{Gray}{wa\_xya}}} }$$ (3.47)

$$\underbrace{ \pi^{m}_{n} }_{ \mbox{{\cmssbx wz\_PiB}}^{\mbox{$\bullet$}} } = \overbrace{ \bigl[ \ln \!\! \underbrace{ p_{s,ijk} }_{ \mbox{{\... ...xy\_PsB}}} } \!\! \bigr]^{m}_{n}. }^{ \mbox{{\cmssbx\textcolor{Gray}{w\_xy}}} }$$ (3.48)

ここで,

$$\underbrace{U_{ijk}}_{ \mbox{{\cmssbx xyz\_UCosLatB}} } = \overbrace{u_{ijk}}^{ \mbox{{\cmssbx\textcolor{blue}{xyz\_UB}}}... ...}_{\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{xy\_CosLat}}} \!\!\!\!\!\!}$$ (3.49)

$$\underbrace{V_{ijk}}_{ \mbox{{\cmssbx xyz\_VCosLatB}} } = \overbrace{v_{ijk}}^{ \mbox{{\cmssbx\textcolor{blue}{xyz\_VB}}}... ...}_{\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{xy\_CosLat}}} \!\!\!\!\!\!}$$ (3.50)

である.

3.6.3 物理過程の時間変化項

外力として与えられた物理過程の時間変化項    $$ をスペクトル値に変換し, 3.5.3節で得られた非重力波項 と足し合わせる. 以下の式の左辺 $\left(\hspace{1em} \right)^{\rm phy}$に対応する変数は 無く, 非重力波項との和が    $^{\mbox{$\bullet$}}$に上書きされる.

$$\begin{align*}\begin{split}\biggl( \DP{\tilde{\zeta}^{m}_{k,n}}{t} \biggr)^{\rm ... ...ox{{\cmssbx\textcolor{Gray}{wa\_DivMu\_xya}}} } \biggr), \end{split}\end{align*}$$ (3.51)

$$\begin{align*}\begin{split}\biggl( \DP{\tilde{D}^{m}_{k,n}}{t} \biggr)^{\rm phy}... ...ox{{\cmssbx\textcolor{Gray}{wa\_DivMu\_xya}}} } \biggr), \end{split}\end{align*}$$ (3.52)

$$\biggl( \DP{\tilde{T}^{m}_{k,n}}{t} \biggr)^{\rm phy} = \overbrace{ \biggl[ \!\! \underbrace{ \biggl( \DP{T_{ijk}}{t} \... ...tPhy}}} } \!\! \biggr]^{m}_{n}, }^{ \mbox{{\cmssbx\textcolor{Gray}{wa\_xya}}} }$$ (3.53)

$$\biggl( \DP{\tilde{q}^{m}_{k,n}}{t} \biggr)^{\rm phy} = \overbrace{ \biggl[ \!\! \underbrace{ \biggl( \DP{q_{ijk}}{t} \... ...tPhy}}} } \!\! \biggr]^{m}_{n}. }^{ \mbox{{\cmssbx\textcolor{Gray}{wa\_xya}}} }$$ (3.54)

ここで,

$$\underbrace{ \biggl( \DP{U_{ijk}}{t} \biggr)^{\rm phy} }_{ \mbox{{\cmssbx xyz\_DUDtPhyCosLat}} } = \overbrace{\biggl( \DP{u_{ijk}}{t} \biggr)^{\rm phy}}^{ \mbox{{... ...}_{\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{xy\_CosLat}}} \!\!\!\!\!\!}$$ (3.55)

$$\underbrace{ \biggl( \DP{V_{ijk}}{t} \biggr)^{\rm phy} }_{ \mbox{{\cmssbx xyz\_DVDtPhyCosLat}} } = \overbrace{\biggl( \DP{u_{ijk}}{t} \biggr)^{\rm phy}}^{ \mbox{{... ...}_{\!\!\!\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{xy\_CosLat}}} \!\!\!\!\!\!}$$ (3.56)

である.

3.7 セミインプリシット時間積分
(サブルーチン TimeIntegration)

注: 本節は書き換え中であるため, 支配方程式系とその離散化の文章とは対応していない.

3.3節 述べた時間積分の2つの手続きは TimeIntegration で行われる. ここでは TimeIntegration 内の計算の流れを, 物理量とモデ ル変数の対応を示しながら説明する.

時間積分の第1段階は (3.1) を解くことである. 時間に依存しない行列 $\tilde{\underline{M}}^{m}_{n}$ はサブルーチン SemiImplMatrix で計算され, LU分解されて    $$として TimeIntegration で使用される. 本来SemiImplMatrixは, サブルーチンDynamicsの冒頭で呼ばれるが, ここに $\tilde{\underline{M}}^{m}_{n}$ を表す式とモデル内の変数との対応を示すこととする.

$$\underbrace{\tilde{\underline{M}}^{m}_{n}}_{ \mbox{{\cmssbx\textcolor{PineGreen}{wzz\_siMtxM}}} } \equiv ( 1-2 \!\!\!\!\!\! \overbrace{\Delta t}^{ \mbox{{\cmssbx\t... ...^{ \mbox{{\cmssbx\textcolor{PineGreen}{wz\_DisCoefM}}} } \!\!\!\!\!\!\!\!\!\! )$$

$$\ \ \ \ - ( \Delta t )^{2} [ \underbrace{ \overbrace{\underline{W... ...\!\!\! \mbox{{\cmssbx\textcolor{PineGreen}{wz\_LaplaEigVal}}} \quad}\! \biggr).$$ (3.57)

この $\tilde{\underline{M}}^{m}_{n}$を用い, $\tilde{\Dvect{f}}^{m}_{n}$ を求める. () を再掲し, 各項の下にモデル内の変数名を記す.

$$\begin{align*}\begin{split}\underbrace{\tilde{\Dvect{f}}^{m}_{n}}_{ \mbox{{\cmss... ...r.2ex\hbox{]}} }_{ \mbox{{\cmssbx w\_siPi}} } \Biggr\} . \end{split}\end{align*}$$ (3.58)

(3.1)は サブルーチン LUSolve により解かれる. 結果は $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$として 変数    $$に格納される. LUSolve の引数は, $\tilde{\Dvect{f}}^{m}_{n}$および LU 分解された $\tilde{\underline{M}}^{m}_{n}$ とそのピボット (    $$,    $$,    $$) である.

3.7.1 $t+\Delta t$の値の算出

(3.5) と (3.6) を解く. 以下に式とモデル変数との対応を記す.

$$\underbrace{\tilde{\pi}^{m,t+\Delta t}_{n}}_{ \mbox{{\cmssbx w\_PiA}}^{\mbox{$\bullet$}} } = \underbrace{\tilde{\pi}^{m,t-\Delta t}_{n}}_{ \mbox{{\cmssbx\te... ...{ \mbox{{\cmssbx wz\_siDivAvrTime}} } \Bigg] }_{ \mbox{{\cmssbx w\_siDPiDt}} },$$ (3.59)

$$\underbrace{\tilde{\Dvect{\zeta}}^{m,t+\Delta t}_{n}}_{ \mbox{{\cmssbx wz\_VorA}}^{\mbox{$\bullet$}} } = \left( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{M} }... ...}}{t}^{\rm NG}}_{ \mbox{{\cmssbx\textcolor{blue}{wz\_DVorDtN}}} } \Biggr\} \, ,$$ (3.60)

$$\underbrace{\tilde{\Dvect{D}}^{m,t+\Delta t}_{n}}_{ \mbox{{\cmssbx wz\_DivA}}^{\mbox{$\bullet$}} } = 2 \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} - \underbrace{\tilde{\Dvect{D}}^{m,t-\Delta t}_{n}}_{ \mbox{{\cmssbx\textcolor{blue}{wz\_DivB}}} },$$ (3.61)

$$\underbrace{\tilde{\Dvect{T}}^{m,t+\Delta t}_{n}}_{ \mbox{{\cmssbx wz\_TempA}}^{\mbox{$\bullet$}} } = \Dinv{1-2\Delta t \tilde{\cal D}_{H,n}^{m}} \Biggl\{ \underbrac... ...ect{D}}^{m}_{n}}^{t} \Biggr] }_{ \mbox{{\cmssbx wz\_siDTempDt}} } \Biggr\} \, ,$$ (3.62)

$$\underbrace{\tilde{\Dvect{q}}^{m,t+\Delta t}_{n}}_{ \mbox{{\cmssbx wz\_QVapA}}^{\mbox{$\bullet$}} } = \Dinv{1-2\Delta t \tilde{\cal D}_{q,n}^{m}} \Biggl\{ \underbrac... ...\! \mbox{{\cmssbx\textcolor{blue}{wz\_DQVapDtN}}} \!\!\!\!\!\!\!\!\!} \Biggr\}.$$ (3.63)