Subsections

3 時間方向の離散化

1 運動方程式と圧力方程式

空間離散化された運動方程式空間離散化された x 方向運動方程式, 空間離散化された z 方向運動方程式と圧力方程式 空間離散化された圧力方程式を時間方向に離散化する. 音波に関連する項は短いタイムステップ で離散化し, その他 の項は長いタイムステップ で離散化する. 音波に関連する項の離 散化には HE-VI 法を採用し, の式は前進差分, の式は後退差分 (クランク・ニコルソン法)で離散化する. その他の項の離散化にはリープフロッ グ法を用いる. 離散化した式の計算はまず の式から行う. 得られた を用いて を計算し, を用いて を計算する.

運動方程式の各項のうち, 音波に関係しない項を として まとめると, 運動方程式と圧力方程式は以下のように書ける.

    (56)
    (57)
    (58)

ただし の式には音波減衰項 を加えてある (Skamarock and Klemp, 1992). 音波に関連しない項 は,
    (59)
    (60)

であり, 時刻 で評価することにする. 但し, 中心差分でリープフロッグ法を用いるため, 数値粘性項 Diff を追加してある.

1 音波に関連する項の時間方向の離散化

1 水平方向の運動方程式の離散化

uwpi:uを時間方向に離散化すると以下のようになる.

    (61)

2 鉛直方向の運動方程式と圧力方程式の離散化

HE-VI 法を用いるので, の式を連立して解く. の式におい て音波減衰項は前進差分, 圧力項は後退差分で離散化する. の式にお いて水平微分項はuwpi:u_sabunで求めた を用いて離散化し, 鉛直微分項は後退差分で離散化する.

    (62)


     
    (63)

ここでは簡単のため格子点位置を表す添字は省略した. uwpi:pi_sabun 式に uwpi:w_sabun を代入して を消去する.
 
   
     
    (64)

uwpi:sabun 式右辺を空間方向に離散化し, 格子点位置を表す添字を付けて表すと以下のようになる (計算の詳細は appendix-a 参照).
     
     
     
     
     
     

但し平均場の量は鉛直方向にしか依存しないので 方向の添字のみ 付けてある.

3 境界条件

上下境界を固定壁とする場合, 境界条件は上部下部境界で,

    (65)
    (66)

である.

下部境界:

下部境界()について考える. この時 uwpi:w_sabun 式に 添字を付けて書き下すと,

 
  (67)

となる. したがって uwpi:sabun_ik 式は以下のようになる.
     
   
     
    (68)

上部境界:

上部境界()について考える. この時 uwpi:w_sabun 式 を添字を付けて書き下すと,

 
  (69)

となる. したがって uwpi:sabun_ik 式は以下のようになる.
     
     
   
     
    (70)

4 圧力方程式の時間積分方法

uwpi:sabun_ik, uwpi:sabun_kabu, uwpi:sabun_joubu 式を連立すると, 以下のような行列式の形式で書く ことができる.

     
    (71)

この連立方程式を解くことで を求める. この連立方程式の係数は以下の ように書ける.
 
     
 
 
 
     
 
     
 
     
 
     
 
     

ただし,
     


 
     

である.

2 音波に関連しない項の時間方向の離散化

運動方程式の音波に関連しない項 uwpi:u, uwpi:w 式を 離散化する.

    (72)
    (73)

ここで, Adv は移流項, D は粘性拡散項, Buoy は浮力項, Diff は数値粘性項である. それぞれの項を書き下すと,
(74)
(75)

であり, 浮力項は,
 
     
    (76)

であり, 粘性拡散項は,
 
     
    (77)
 
     
    (78)

である. 数値粘性項は,
(79)
(80)

である. は乱流エネルギーの時間発展方程式から計算し(詳細は後述), は以下のように定める.
    (81)
    (82)

ここで は水平・鉛直方向の格子間隔を意味し, はそれぞれ,
    (83)

とする.

2 熱力学の式と混合比の保存式の離散化

熱の式と混合比の保存式の右辺をまとめて で表し, 時間方向にリープフロッグ法を用いて離散化する.

(84)
(85)
(86)
(87)

ここで,
 
    (88)
 
    (89)
 
    (90)
 
    (91)

である. 移流を中心差分で安定して解くために, 数値粘性項 Diff を追加してあ る. また, 項は湿潤飽和調節法より決めるため, それらの項を含めない.

, , , をまとめて で表し, それぞれの項を書き下す. 移流項は,

(92)

であり, 基本場の移流項は,
    (93)

である. 粘性拡散項は CReSS と同様に 1.5 次のクロージャーを用いることで,
 
    (94)

となり, 基本場の粘性拡散項は,
(95)

となる. 数値粘性項は,
(96)

である. は乱流エネルギーの時間発展方程式から計算する(詳細は後述). は nu 式を利用する.

凝縮加熱項

(97)

である.

散逸加熱項

(98)

と与える. ここで である.

放射強制 は計算設定ごとに与える.

雲水から雨水への変換を表す , は以下のようになる.

    (99)
    (100)

雨水の蒸発を表す は以下のようになる.
    (101)

降水による雨水フラックスを表す は以下のように書ける.
    (102)
    (103)

1 湿潤飽和調節法

Klemp and Wilhelmson (1983), CReSS ユーザーマニュアル(坪木と榊原, 2001) では, 水蒸気と雲水の間の変換を表す は, Soong and Ogura (1973) において開発された 湿潤飽和調節法を用いる. この方法は の断熱線と, の 平衡条件( は化学ポテンシャル)の交わる温度・圧力・組成を 反復的に求める数値解法である. 以下ではそのやり方を解説する.

1 飽和蒸気圧を用いる場合

湿潤飽和調節法を用いる場合, まず始めに risan:time-div_theta - risan:time-div_qr 式から求まる量に を添付し, , , , とする. 水に対する過飽和混合比

    (104)

, もしくは雲粒混合比が なら ば, 次式を用いて暫定的に , , を求める.
(105)
(106)
(107)

ただし, である. もしも ならば, 暫定的に得られた値を 付き のものに置き換え, moistajst_theta1 - moistajst_qc1 式 の値が収束するまで繰り返し適用する. 普通, 高々数回繰り返せば収束し, 調整後の値が得られるそうである.

もしも の場合には,

    (108)
    (109)
    (110)

とし, 繰り返しを中止する.

2 圧平衡定数を用いる場合

硫化アンモニウムの生成反応

    (111)

のような, 2 種類の気体 1 モルづつから凝縮物質 1 モルが 生成されるような生成反応の場合の, 湿潤飽和調節法を考える.

硫化アンモニウムの生成反応の圧平衡定数は,

    (112)

である. 圧平衡定数を用いることで, 任意の温度に対する アンモニアと硫化水素のモル比の積を求めることができる.

任意の温度 における NHSH の生成量を とすると, 圧平衡定数の式は以下のように書ける.

     
    (113)

解の公式を使うと, 生成量 X は以下となる.
     
    (114)

根号の符号は の場合にとりうる の値を 仮定することで決める. の場合, 明らかに
    (115)

である. ここで木星大気を想定し, であることを仮定すると である. そこで def_X_NH4SHの根号の符号は のとき となるよう, 負を選択する.
    (116)

の満たすべき条件は,
    (117)

である. 上記の条件を満たさない場合には とする.

が NH4SH-condition 式の条件を満たすならば, 次式を用いて暫定的に , , を求める.

    (118)
    (119)
    (120)
    (121)

ただし, であり, はそれぞれ, 生成量 に対応する NH と HS の混合比である. 温位が収束するまで反復改良を行う.

3 乱流運動エネルギーの式

Klemp and Wilhelmson (1978) および CReSS (坪木と榊原篤志, 2001) と同様 に, 1.5 次のクロージャーを用いる. 乱流エネルギーの時間発展方程式 をリープフロッグ法を用いて時間方向に離散化すると, 以下のようになる.

    (122)

ここで,
 
    (123)

である. CReSS にならい, 移流項を で, 移流項以外を で評価した.

に含まれる各項は以下のように書き下すことができる.

(124)
(125)
 
     
    (126)
 
    (127)
(128)

ここで , 混合距離 とする. また は以下で与えられる.
(129)
(130)

ただし,
(131)

である.

4 時間フィルター

リープフロッグ法を用いたことによって生じる計算モードの増幅を抑制するた め, Asselin (1972) の時間フィルターを長い時間刻みで 1 ステップ計算する 毎に(実際には短い時間刻みの計算を ステップ計算する毎に)適用する.

たとえばuwpi:u_sabunを用いて を計算する場合, 以下のように時間フィルターを適用する.

 
     
(132)

ここで はフィルターの係数であり, その値は 0.05 を用い る. uwpi:w_sabun, uwpi:pi_sabunの計算に対しても同様 に時間フィルターを適用する.

5 スポンジ層

境界面付近での波の反射を抑えるために, 基礎方程式の付加的な項を付け加える.

    (133)

ただし, は任意の予報変数であり, は客観解析値等の既知の 値である. この項は1 つ前のタイムステップ で計算され, 小さいタイムステップで扱われる予報変数に対しても, 移流項や数値粘性項と同様に の大きなタイムステップ間の値とし て評価される。具体的には,
    (134)
    (135)
    (136)
    (137)

とする. 但し はエクスナー関数の基本場である.

はそれぞれ水平方向には各境界面に向かって, 鉛直 方向には上境界面に向かって小さくなる減衰係数である. これらの減衰係数は, 水平方向には吸収層の厚みを とし, の範囲を とすれば,

     
     
    (138)

であり, 鉛直方向には吸収層の厚さを とし, の範囲を とすれば,
     
    (139)

である. ここで, はそれぞれ水平・鉛直方向の減衰定数 である. は時間の逆数の次元を持ち, それらの逆数 は e-folding time と呼ばれる. e-folding time は通常 100 - 300 s に設定する. また吸収層の厚み はそれぞれ, 水平方向には数格子分, 鉛直方向には上面から1/3 程度設定すれば良い.

SUGIYAMA Ko-ichiro 2011-06-14