数値予報モデルの説明

Return Home  短期予報   １週間～１ヶ月予報  　
数値予報モデルの説明　操作のヒント　検討事項（CFL条件、誤差）　問題点・改善点 
Lorenz Strange Attractor  初期値鋭敏性の検討図  誤差の増大図

Last updated 2002-1/14,1/16,1/17,1/26,1/30,1/31,2/2,7/13,7/14,7/15,2006-3/05

本HPにおける数値予報プログラムは、気象学における「プリミティブ方程式」の数値解をグラフ化
して見るものです。　値を適宜キーインし「Draw」ボタンを押すと、グラフが描かれます。

以下、この数値予報モデルの説明を行います。また、検討事項は、随時 update します。

この数値予報プログラムは、下記プリミティブ方程式：

運動方程式（大気塊１ｋｇを質点として考え、力の釣合を表す）：
　　①∂ｕ／∂ｔ＝－（ｕ∂ｕ／∂ｘ＋ｖ∂ｕ／∂ｙ＋ｗ∂ｕ／∂ｚ）　
　  　　　　　　　＋２Ω(ｓｉｎφ)ｖ－(1／ρ)∂ｐ／∂ｘ＋Ｆx　
　　②∂ｖ／∂ｔ＝－（ｕ∂ｖ／∂ｘ＋ｖ∂ｖ／∂ｙ＋ｗ∂ｖ／∂ｚ）　
  　　　　　　　　－２Ω(ｓｉｎφ)ｕ－(1／ρ)∂ｐ／∂ｙ＋Ｆy　
　　③０＝－∂ｐ／ρ∂z－ｇ　　
連続の式（大気１ｋｇの質量保存則、密度変化を表す）：
　　④∂ρ／∂ｔ＝－（ｕ∂ρ／∂ｘ＋ｖ∂ρ／∂ｙ＋ｗ∂ρ／∂ｚ）
  　　　　　　　　－ρ（∂ｕ／∂ｘ＋∂ｖ／∂ｙ＋∂ｗ／∂ｚ）
熱力学の第一法則（大気塊１ｋｇ中の温位の増分）：
　　⑤∂θ／∂ｔ＝－（ｕ∂θ／∂ｘ＋ｖ∂θ／∂ｙ＋ｗ∂θ／∂ｚ）＋Ｈ
水蒸気量保存の法則（大気塊１ｋｇ中の水蒸気量の増分）：
　　⑥∂ｑ／∂ｔ＝－（ｕ∂ｑ／∂ｘ＋ｖ∂ｑ／∂ｙ＋ｗ∂ｑ／∂ｚ）＋Ｍ
気体の状態方程式：
　　⑦Ｐ＝ρＲＴ

のうち、①②の風に関して、次のような漸化式（オイラー法）に置き換え、モデルとしたものです。
（その他の式は一切無視しました。今後組み込む予定です。）
（なお、時間ができれば、計算精度向上のため、ルンゲ・クッタ法もtryするつもりです。）

xn1=xn+(-(xn*gradxnx+yn*gradxny)+2*omega*Math.sin(phi)*yn-(1/dens)*(gradpresx)+fx)*dt;
yn1=yn+(-(xn*gradynx+yn*gradyny)-2*omega*Math.sin(phi)*xn-(1/dens)*(gradpresy)+fy)*dt;

上式に下記の値を代入し、繰返し計算し、所要時刻の予想値を算出します。
速度の傾斜量や気圧の傾斜量は気象の変化に密接に関係する量であり、計算により算出
すべきですが、ここでは、初めての試みなので、簡略化し一定値としました。但し、
それらの値は手入力にて変更可能です。なお、予め与えた値は、必ずしも気象現象と
して起こり得るものとは限りません。そのためにも、起こり得るであろう値を様々に
与えてシミュレートすると良いでしょう。
------------------------------        -----------------------------------
変数名　：　意味　　　　　　　　　　　　　スクリーン上の表示と初期値  （単位等)
-------：--------------------        -----------------------------------
xn,xn1 ：   風速の東西成分ｕ             　u0=　　　　　20.0 　　(m/s) 
yn,yn1 ：   風速の南北成分ｖ             　v0=　　　　　10.0 　　(m/s) 

gradxnx：　 uのｘ方向の変化率     　       ux=      0.00001　   (m/sm) 
gradynx：   ｖのx方向の変化率        　    vx=      0.00001   　(m/sm)
　　　　　　　　　　　　　　　　　　　　　(1000kmあたり10m/sの傾斜としました)
       
gradxny：　 uのy方向の変化率               uy=     0.00001  　 (m/sm)        
gradyny：   vのy方向の変化率               vy=     0.00001   　(m/sm)  
　　　　　　　　　　　　　　　　　　　　　(1000kmあたり10m/sの傾斜としました)

gradpresx： 気圧のx方向の変化率         　　px=      0.001 　   (Pa/m) 
gradpresy： 気圧のy方向の変化率      　     py=      0.001 　   (Pa/m)
　　　　　　　　　　　　　　　　　　　　　(1000kmあたり10hPaの傾斜としました)

fx，fy：    ｘ方向、ｙ方向の摩擦力     摩擦力FxとFy=-0.001(地表面付近の値は大きい)。
　　　　　　　　　　　　　　　　　　　　摩擦力は、加速度方向とは逆向きと考えられるため、
　　　　　　　　　　　　　　　　　　　　マイナス表示にしました。しかしどれ位のオーダー
　　　　　　　　　　　　　　　　　　　　の数字か分かりません。なお、自由大気の高度に
　　　　　　　　　　　　　　　　　　　　おいては、０と考えたいですが、必ずしもそうでは
　　　　　　　　　　　　　　　　　　　　ない様です。

                                        なお、２００２．７／１３の神奈川気象予報士会でこの数値予報ＰＧＭ
　　　　　　　　　　　　　　　　　　　　を説明した際、同会参加者の一人から、摩擦力は、速度ベクトルと逆向きに
　　　　　　　　　　　　　　　　　　　　作用するので、摩擦項の符号は速度に対応して正負を設定すべき、との
　　　　　　　　　　　　　　　　　　　　有り難いご指摘をいただきました。近日中にプログラムを修正予定です。　 



fx，fyの値は、　［　気圧傾度力Ｐ＝コリオリ力Ｃ＋摩擦力Ｆ　］の力の釣り合いから求めるものとしよう。

Cをコリオリ力、等圧線と風速のなす角度をαとして、摩擦力Ｆ=コリオリ力Ｃ＊tanα　
ここにおいて、コリオリ力Ｃを速度の成分方向(x,y)に分解すればよいだろう。
すなわちコリオリ力を、(u,v)方向に分解しよう。コリオリパラメーターをｆ(=2Ωsinφ)と置く。
Ｃのu方向Cxは、Cx=fu
Ｃのv方向Cyは、Cy=fv
となる。
すなわち、
fx=futanα=(2Ωsinφ)u(tanα)
fy=fvtanα=(2Ωsinφ)v(tanα)
となる。

以上の式で、fx,fyを計算することとしよう。
なお、摩擦力の符号は速度の符合と反対に設定します。
プログラムでは、次の式となる：
fx=(2*omega*sin(phi)*(-xn)*(tanα))
fy=(-2*omega*sin(phi)*(-yn)*(tanα))


なお、tanαは、気圧及び風速のx,y成分を考えることにより、
等圧線とｘ軸との角度をβ、速度Ｖとｘ軸との角度をγとするとき、
α＝γ－β
     sinγ=v/Ｖ　
     cosγ=u/Ｖ　

     sinβ=px/Ｐ
     cosβ=py/Ｐ
となる関係を使って、
tanα=sinα/cosα
     =sin(γ－β)/cos(γ－β)
     =(sinγcosβ－cosγsinβ)/(cosγcosβ+sinγsinβ)
     =(v*py-u*px)/(u*py+v*px)
と表される。このtanαを上のfx,fyに代入します。
これで、摩擦力に関係した部分の数値計算の準備は整ったと考えられます。

以上の  摩擦力を考慮したプログラム を作成中です。検討のために、時間がかかります。

omega：  　地球の自転の角速度　　　　　Ω=0.000073　(rad/s)　をプログラム内で与えています。　
phi：　　　 予想対象地点の緯度φ　　　　緯度phi=35　　(度）（プログラム内でradへ変換します）
dens：　　  大気の密度                 dens=1.05（kg/m3）(自由大気となる850hPa付近の値に固定しました)

dt：        積分時間間隔dt分               dt=      6.0　　(min)
            計算１回あたり現象の発達・進行する時間 　　
ｎ:         計算打切り回数n                n=　　　　360　　　　
            dt＊n　により予想対象時刻を指定します。例えば、
            6分＊360回→2160分＝36時間先までの予想値を求めます。　
--------------------------------------------------------------------------------------------
検討事項を参照してください。
AZURE トップページへ戻る Hajime satoh,  All Rights Reserved. 2002-2006
ＡＺＵＲＥ 青空