気象基礎演習Ⅱ
可視化ツールを使った解析データの利用と、気象を支配する物理法則の体験的理解

解析データとは解析データとは何らかの合理的な方法に基づいて推定された時間的・空間的に一様なデータセット。
特に数値モデルを使って物理法則から推定したデータは同化データと呼ばれる。
また、解析データや数値予報の予報値は空間的に一様な格子点上に並んでいるのでGPV（ Grid
Point Values ）とも呼ばれる。

左図は2008年6月1日00Zに観測された高層データ（TEMP報）の分布（赤丸）を表しているが、空間的に著しい偏りが見られ、このままでは定量的な分析に不向きである。
図の中の格子は緯度・経度共に2.5度で引かれているが、これは代表的な全球の解析データの格子間隔にほぼ等しい。つまり、解析データでは図の網の目のように、地球全体を海陸分布に関係なく一様に覆っているため、様々な定量的分析に使用できる。
2008年6月1日00Zに観測された高層データ（TEMP）の分布（赤丸）解析データ（GPV）の利点：量的分析が容易
dV

1=− f k ×V − ∇ pdt

大気の運動方程式は

この中の気圧傾度ベクトル ∇ p=



∂p
∂x

,

∂p
∂y



解析データでは空間差分で大雑把に見積もれる。

は

例えば、左下図の原点Oにおける気圧傾度は


∂p

観測データ

y

∂x O

解析データ

P4

PA

PD

PB
P2

P1

≈

P B −P A
2 x

,


∂p

∂y O

≈

P D −P B
2 yから近似的に求められる。
一方、空間的に不均一に分布した観測データ

P3

O

V ：風速ベクトル水平成分
f ：コリオリパラメータk ：鉛直方向の単位ベクトルr：密度、 p ：気圧と表現できる。

P 1 , P 2, P 3, P 4xでは見積もりが困難。

PC

y

x

 
∂p

∂x O

≈? × P 1? ? ×P 2? ? ? ×P 3? ? ? ? ×P 4

GPV可視化ソフトの必要性
人間コンピュータ

演算能力

低い

非常に高い

画像認識能力高い

低い

解析データ（GPV）はコンピュータを使った数値的なデータ処理には非常に便利だが、そのままでは大量の数値データの集まりであり、そのままでは人間に理解不能。

GPV可視化ソフトの必要性
コンピュータの高い演算能力を利用して、数値情報→図情報に変換。
人間の画像認識能力は高いので、大量の数値情報も容易に把握できる。

利用目的に応じて様々なGPV可視化ソフトが作られている。
●

１次元データ、２次元データ、３次元データ、４次元データ●

静止画像、動画●

ラスターグラフィックス、ベクトルグラフィックス

様々なGPV可視化ソフト

低レベル

OSや特定のプログラム言語に依存
GTK, Win32API, Java VM 。線や多角形・文字等を描く単純な機能しかない。

中レベル

等値線や陰影など高級なライブラリ関数はあるものの、特定のプログラ
OpenGL, PGPLOT, NCARG
ム言語に依存し、ソースプログラムの作成とコンパイルが必要。

高レベル

独自のコマンド群で構成され、特定のプログラム言語を必要としない。
利用形態は対話型でのコマンド入
力か、スクリプトによるバッチ処理。

gnuplot, grads, GMT

gnuplot, GrADS, GMTの比較
gnuplot

GrADS

GMT

コマンド逐次入力によるコマンド逐次入力によるインターフェース 対話型
対話型
スクリプトによるバッチ処理
スクリプトによるバッチ処理 スクリプトによるバッチ処理
描画機能
入力データ

やや低い

中程度

高機能

テキスト、バイナリ

バイナリ（GrADS標準形式 テキスト、、GRIB、netCDF、BUFR、 バイナリ（netCDF、
HDF、GRIB2、Station） その他複数の形式）

地図（海岸線） なし

３段階の解像度の地図デ ５段階の解像度の地図デ
ータを持つ
ータを持つ

地図の投影法 なし

等緯度経度、モルワイデ、 等緯度経度、モルワイデ、
ポーラーステレオ、
ポーラーステレオ、
ランベルト等各種
ランベルト等各種

画面（X11等）、GIF、JPEG、
出力画像の形式
PNG、TeXDraw、
Postscript等各種画面（X11等）、GIF、PNG、Postscript

Postscriptのみ

アニメーション 不可

可能

不可

gnuplot, GrADS, GMTのサンプル画像

gnuplot
GMT

GrADS

GrADSとは● Grid Analysis and Display Systemの略。
● ４次元空間（東西、南北、上下、時間）のデータを自由に切り取って表示させる。
● 様々なOS（Windows、MacOSX、Linux、FreeBSD、商用UNIX）で利用可能。

自分のPC等にインストールする方法● 操作はコマンドを入力し画面表示を見ながら対話型で行う。
● 制御文・関数を備えた独自のスクリプトでバッチ処理も可能。
● 作成された画像はファイルに出力できるだけでなく、端末に動画出力も可能。
● データはGPVのみを保存したデータファイルと、データの時空間構造やデータの種類を記述したコントロールファイルの2つから構成される。
● 標準形式のデータファイルはC言語やFORTRANなどのプログラム言語

で直接読み書きすることが容易。
● データファイルは標準形式以外に、GRIB、GRIB2やnetCDF、BUFRなど

気象で使われているファイル形式にも対応。

GrADSコントロールファイル/home/kisoen2/grads/sample_data/jra1st10dJAN08.ctl

GrADSコントロールファイルの表示

GrADSデータの時空間構造と物理量の種類・記録された順序を記述

物
理
量
の
名
前
、
コ
マ
ン
ド
で
表
示
さ
れ
る
と
き
に
使
わ
れ
る
。

DSET
UNDEF

^jra1st10dJAN08.gdt
未定義値-1.0E30

XDEF
YDEF

144 LINEAR
0.0 2.5
73 LINEAR -90.0 2.5

GrADSデータファイルの名前「^」はデータファイルが同じ
ディレクトリにあることを示す。
経度方向144点、東経0度から2.5度間隔
緯度方向73点、南緯90度から2.5度間隔

ZDEF
17 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70
鉛直17層：1000～10hPa：「LEVELS」は不等間隔を表す。
50 30 20 10
TDEF
OPTIONS

20 LINEAR 00Z01JAN2008 12hr
little_endian yrev

VARS
gph

17

6
99

geopotential height (gpm)

tmp

17

99

air temperature (K)

sph
u

12
17

99
99

specific humidity (kg/kg)
zonal wind (m/s)

v
psl

17
0

99
99

meridional wind (m/s)
sea level pressure (Pa)

時間は0Z1JAN2008から12時間で20個「yrev」はデータの南北方向の並びを逆にする意味（つまり北から南へ）

物理量は6種類

ENDVARS
コントロールファイルの終わりを表す、必須

鉛直17層
鉛直12層
コメント
2次元物理量時間の単位
mn：分、hr：時間
dy：日、mo：月、yr：年GrADSデータファイル/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt（例:jar1st10dJAN08.gdt）

データの形式：GPVが4バイト浮動小数点形式で単純に並んでいる。
データの並び：（（（（（GPV、東→西）、南→北）、下→上）、変数の種類）、時刻）
１レコードサイズ=4*144*73=42048バイト

ファイルの先頭gph(1000hPa, 00Z01JAN08)u(10hPa, 00Z01JAN08)v(1000hPa, 00Z01JAN08)gph(10hPa, 00Z01JAN08)
tmp(1000hPa, 00Z01JAN08)v(1000hPa, 00Z01JAN08)psl(00Z01JAN08)

tmp(10hPa, 00Z01JAN08)

最初のレコード
gph(0E,90N)
gph(2.5E,90N)

gph(1000hPa, 12Z01JAN08)gph(2.5W,90N)
gph(0E,87.5N)
gph(2.5E,87.5N)

gph(10hPa, 12Z01JAN08)gph(0E,90S)

sph(1000hPa, 00Z01JAN08)

sph(10hPa, 00Z01JAN08)
u(1000hPa, 00Z01JAN08)

繰り返し続くgph(2.5W,90S)

GrADS利用の実際
1. 端末からGrADSコマンドを入力する
端末で「grads」とコマンドを入力すると、「Landscape mode? ('n' for portrait):」というプロンプトが出るので、横長画面であれば単純にリターンキーを縦長画面であれば「n」を入力する。

2. GrADSファイルを開く（openコマンド）
GrADSからコマンド入力を促すプロンプト「ga->」が出るのを待つ。プロンプトが出たら、目的のGrADSコントロールファイルが「/ABC/defg/hijk.ctl」のとき、コマンド
open /ABC/defg/hijk.ctl
を入力する。コントロールファイルが正しく記述され、データファイルと齟齬がなければ、コントロール
ファイルとデータファイルの名前、データの四次元構造（経度、緯度、高度、時刻）が表示された後、プロンプトが出てコマンド入力を待つ（コントロールファイルやデータファイルが存在しないか誤りが見つかった場合はエラーメッセージが表示される）。

3. 画像表示するための各種パラメータを設定する（setコマンド）
GrADSでGPVを表示するために、地図の投影法や、四次元空間から切り出す断面のパラメータ設定などを「setコマンド」を入力して行う。
set <パラメータのリスト>

4. 画像表示する物理量を指定し画面に表示させる（dコマンド）表示させる物理量の名前が「var123」（コントロールファイルのVARS欄以下に記述された名前）であったときは、コマンド
d var123
で画像表示できる。画像表示出来るものは物理量自身以外に、物理量の和・差・積・商や幾つかの関数で変換した量などがある。

GrADS利用の実際（続き）
5. 作成した画像をファイルに保存する（gxprintコマンド）
GrADSで描画した画像は「gxprint」コマンドで画像ファイルとして保存できる。
ベクトル画像ファイルであるEncapsulated Postscript (eps) ファイルに変換する場合は、コマンドの後にファイル名を指定してgxprint （ファイル名）.epsとする。
一方、GIFやPNG等のラスター画像ファイルに保存する場合は更に解像度を指定してgxprint （ファイル名）.png x1024 y768
等とする。
Ver 2.0 以前での方法

6. GrADSを終了する（quitコマンド）
一連の作業を終えてGrADSを終了させたいときはコマンド
quitを入力するだけでよい。

GrADSの各種設定（setコマンド）●次元の設定GrADSの次元は格子点の番号（X, Y, Z, T）或いは実際の経度・緯度・高度・時刻（LON, LAT, LEV, TIME）の2種類のパラメータで指定できる。
例えば、jra1st10dJAN08.ctlの場合、
X=57, Y=23, Z=6, T=5 は LON=140E, LAT=35N, LEV=500hPa, TIME=0Z03JAN2008に対応する。
東経120度～160度、北緯20度～50度、500hPa、12Z05JAN2008の部分を設定する場合はset lon 120 160
set lat 20 50
set lev 500
set time 12z5jan2008という一連のコマンドを入力すればよい。

●投影法の設定
GPVを地図と共に水平二次元平面に表示させる場合は何らかの投影法を設定する必要がある（球面の地球⇒平面図）。その設定には投影法パラメータmprojを設定する。
set mproj latlon
等緯度経度投影（デフォルト）
set mproj nps
ポーラステレオ投影（北半球側）
set mproj sps
ポーラステレオ投影（南半球側）
set mproj lambert
ランベルト投影set mproj mollweide
モルワイデ投影set mproj off
地図・投影ともになし

GrADSの各種設定（setコマンド）●グラフィックタイプの設定
GrADSの四次元データでグラフィック表示できるのは、その中から「set」コマンドで切り出した二次元以下の部分である。グラフィックタイプはパラメータgxoutを設定して指定する。
▼二次元データ
set gxout contour 等値線による表現（二次元のデフォルト）
set gxout shaded 塗りつぶしによる表現
set gxout vector 矢印による表現（データは東西・南北の2成分が必要）
set gxout barb
矢羽根による表現（データは東西・南北の2成分が必要）
set gxout stream 流線による表現（データは東西・南北の2成分が必要）
set gxout fgrid
値により色分けした格子内の塗りつぶしによる表現
set gxout scatter 散布図による表現（2要素のデータが必要）▼一次元データ
set gxout line
折れ線による表現（一次元のデフォルト）
set gxout bar
棒グラフによる表現
set gxout errbar 誤差幅の表現※三次元データのうち、一つが時間であれば、アニメーションで表示される。
（空間１次元の場合のアニメーション機能のon, off

set looping on, off）

●グラフィックオプションの設定
各グラフィックタイプ毎に細かな設定をするオプションがあり、これも「set」コマンドで設定できる。
例：「gxout」が「contour」の場合、等値線の間隔は「set cint」で、等値線の色は「set ccolor」で、等値線の太さは「set cthick」で、等値線のラベルを付けるかどうかは「set clab」で、それぞれ設定できる。
グラフィックオプションの設定は gadoc/advdisplay.html を参照せよ。

データの表示（dコマンド）「set」コマンドで各種設定を終了後は、「d」に続き表示する物理量の名前（コントロールファイルの「VARS」以下の行にある名前）をコマンド入力するだけで、物理量は画面に表示される。例えば、物理量を「var123」とすると

d var123と入力する。但し、「gxout」が「vector」、「stream」、「scatter」に設定されている場合は、2種類の物理量を以下のように同時に指定しなければならない。

d u123;v123
これらの場合、表示される物理量の次元やデータ構造は、「set」コマンドの設定に従う。
例えば、

set lon 120 150
set lat 30 40
set lev 500
set time 0z1jan2008と設定しておくと、「d var123」で画面に表示されるのは、東経120度～150度、北緯30度～40度、高度500hPa、
時刻0z1jan2008での「var123」の2次元分布である。一方、

set lon 140
set lat 35
set lev 1000 100
set time 0z1jan2008 0z5jan2008と設定した場合は、東経140度、北緯35度、高度1000～100hPa、時刻0z1jan2008～0z5jan2008までの時間－高度断面が表示される。

データの表示（dコマンド）（続き）●「d」コマンドでの次元設定「d」コマンドを入力する際に、物理量の名前のみならず次元の設定も同時に出来る。例えば、

d var123(lon=140,lat=35,lev=500)
はset lon 140
set lat 35
set lev 500
d var123と同じ効果である。

●「d」コマンド中での演算「d」コマンドのパラメータとしては、物理量だけでなく、以下のような物理量間の四則演算等も可能である。

d 1004*tmp+9.8*gph

気温と定圧比熱の積に重力加速度と高度の積の和（Dry Static Energy）

d gph(lev=500)-gph(lev=850)
d ave(u,lon=90,lon=180)

500hPaの高度と850hPaの高度の差（層厚）ave関数を使った東経90度から東経180度までの平均値
なお、演算する物理量間の次元構造は完全に一致していなければならない。

GrADSスクリプト
GrADSにはスクリプト言語があり、一連のコマンド入力を自動的に実行してくれるだけでなく、代入文や制御文等独自の文法を備えていて、きめ細やかな描画の設定が可能になる。
0. コマンドの逐次処理
最も基本的な機能で、手で入力するコマンドを、処理する順序に並べておけばよい。但し、コマンドやパラ
メータは前後を引用符「'」で囲んでおくこと（引用符で囲まれた部分はそのまま文字としてGrADSが解釈するが、囲まれていないものは制御文等の予約語もしくは変数と解釈される）。
例えば、手入力でset gxout shaded
set time 0z1jan2008と入力されるコマンドは、スクリプトでは'set gxout shaded''set time 0z1jan2008'というように記述される。
1. 変数と配列
スクリプトでは変数や配列が利用できるので、長いファイル名や数多くのパラメータ入力の繰り返しも簡潔に記述することができる。例えば、GrADSで時刻を次々と変更する場合、コマンド入力でset time 1JAN2008……………………
set time 1FEB2008……………………
set time 1MAR2008

GrADSスクリプト（続き）というように年月日を変えて様々に処理する必要がある時、スクリプトでも同様の記述をすればよいが、破線の部分が全く同じ処理をする場合、以下のように変数・配列を使ってより簡潔に記述できる。

monnam.1='JAN'; monnam.2='FEB'; monnam.3='MAR'
day=1; mon=1; yr=2008
while(mon<4)tm=day % monnam.mon % yr
'set time 'tm…………………
mon=mon+1
endwhile
ここで、青色は変数・配列、黒色は演算子、緑色は制御文の予約語を表す。変数mon=2のとき、配列
monnam.monは「FEB」の値をとるので、変数tmは結局「1FEB2008」が代入され（式の評価の際に変数は代入された文字として展開される。演算子「%」は文字列の結合を意味する。）、'set time 'tm
は'set time 1FEB2008'と解釈される。

GrADSスクリプト（続き）
3. 演算子
変数・配列の演算のために演算子が用意されている。一般的な四則演算用の「+」、「-」、「*」、「/」演算子以外に剰余計算の「%」演算子がある。この場合、演算は全て浮動小数点数として行われる。論理演
算子としては論理和「｜」、論理積「&」、否定「!」等の演算子が用意されている。この他、「>」、「>=」、「<」、「<=」のような比較演算子も使える。
4. 制御文スクリプト実行の流れを制御する制御文としては、条件分岐のための「if」文と、繰り返し実行のための「while」文がある。「if」文の構造はif ( 条件式 )
実行文群
else
実行文群
endif
となっている。一方、「while」文の構造は

while ( 条件式 )実行文群
endwhileとなっている。「while」文では条件式が成立している間、何回でも実行を繰り返すが、「break」文や「continue」文で処理の流れを変えることが出来る。「break」文に達すると、即座に繰り返しのループから抜け出すように制御される。一方、「continue」文に達すると、それ以降の処理を中止し、ループの最初に戻るよう制御される。

GrADSスクリプト（続き）●スクリプトの実例
スクリプトを使うと描画で細かい設定が可能になる。下の図は右にあるスクリプトを使って描いた図で、等値線ごとに色や太さ、ラベルなどを変えている。

z500描画GrADSスクリプト'open /home/kisoen2/grads/sample_data/jra1st10dJAN08.ctl'
zmin=4800; zmax=6000; zint=50; lblint=4
tm='00Z01JAN2008'
等値線の最小値4800、最大値6000、'set time 'tm
間隔50、主コンターの間隔4に設定'set lat 20 90''set lev 500'
'set gxout contour''set mproj nps'
等値線の値をzminに、カウンターを0に設定'set frame off'
'set map 7 1 1'
等値線の値がzmaxになるまで繰り返す
zval=zmin; i=0
while(zval<=zmax)'set clevs 'zval
if(i=0)'set ccolor 5'
カウンターが0のときは水色で'set cthick 6'
太い等値線にラベルをつける'set clab on'
else'set ccolor 1'
カウンターが0以外のときは白色で'set cthick 1'
細いラベルのない等値線'set clab off'
endif
'd gph'zval=zval+zint
等値線の値とカウンターの更新
i=i+1
if(i=lblint)i=0「draw title」の後に続くテキストを
endif
図のタイトルとして適当に表示する。
endwhile
このとき変数「tm」は文字列に展開'draw title Z500 ('tm')'
される。

GrADSデータの読み書き「GrADSデータファイル」、「GrADSコントロールファイル」のところで述べたように、GrADSデータファイルは4バイト不動小数点数（バイナリ）を並べただけの単純な構造をしているため、FORTRANやCなどのプログ
ラミング言語で読み書きすることは比較的容易である。
例えば、GrADSデータファイルjra1st10dJAN08.gdtの00Z05JAN2008における500hPa高度を読む
FORTRANプログラムおよびCプログラムは以下のようなものになる。

CプログラムFORTRANプログラム
program z500_read
implicit none
integer,parameter :: idim=144,jdim=73
integer :: ios,irec
real :: z500(idim,jdim)
open(10,file='jra1st10dJAN08.gdt',iostat=ios, && action='read',recl=4*idim*jdim, && access='direct',form='unformatted')if(ios /= 0) stop 99
irec=4*2*(4*17+12+1)+6
read(10,rec=irec) z500
close(10)end program z500_read

gfortranの場合で、ifortの場合
は4バイトを1単位としているので idim*jdim とする。

#include <stdio.h>#define IDIM
#define JDIM

144
73

int main()
{
int
irec;
float
z500[JDIM][IDIM];
FILE
*fp;
if( ( fp = fopen( “jra1st10dJAN08.gdt”,“r”) )== NULL ) return 99;
irec=4*2*(4*17+12+1)+5;
fseek( fp, IDIM*JDIM*irec*sizeof(float), SEEK_SET );fread( z500[0], sizeof(float), IDIM*JDIM, fp );fclose( fp );return 0;}

GrADSファイルの作成（温位データ）
GrADSは物理量間の簡単な四則演算をする機能を備えているが、複雑な関係式から値は計算できないので、この場合はFORTRANやC等のプログラミング言語で計算して必要なGrADSファイルを作成する必要がある。
この例として、気温と気圧から温位を求め、そのGrADSファイルを作成して、画像に出力してみる。

温位とは温位（Potential Temperature）は気温T 、気圧p の大気では

R





≡T p0/ pCpという式で定義される。ここで、Cp (=1004JK-1kg-1)は乾燥大気の定圧比熱、R (=287JK-1kg-1)は乾燥大
気の気体定数を表す。また、p0 =1000hPaである。
温位は断熱状態にある空気では気塊固有の値を持ち変化しない。非断熱加熱率Q で空気が加熱された
場合、温位は下記の式のように変化する。

d
dt

=

Q
C pT

温位は乾燥大気の静的安定度の目安にもなっていて、ddz0  安定 ,と関係づけられている。

ddz=0  中立 ,ddz

0  不安定GrADSファイルの作成（温位データ）(続き)
GrADSデータ/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt から指定した日付の温位を計算
しGrADSデータファイルを作成する。
FORTRANプログラム
program ptcalc
implicit noneopen(infile,file=finame,iostat=ios,action='read', && recl=lrec,access='direct',form='unformatted')if(ios /= 0) stop 33

integer,parameter :: idim=144,jdim=73,ldim=17
integer :: infile=11,ioutfl=22
integer :: idini(4)=(/2008,1,1,0/),tmax=20
integer :: idout(4)=(/2008,1,5,0/)
integer :: i,j,l,ios,lrec,ntime,irec
real,dimension(ldim) :: plev= && (/1000,925,850,700,600,500,400,300, &&250,200,150,100,70,50,30,20,10/)
real :: cp=1004.0,rgas=287.0,p0=1000.0,kappa,factor
real,dimension(idim,jdim,ldim) :: tmp,pt
character(len=100) :: finame,foname

irec=ntime*(4*ldim+12+1)+ldim
do l=1,ldim
read(infile,rec=irec+l)((tmp(i,j,l),i=1,idim),j=1,jdim)
factor=(p0/plev(l))**kappa
do j=1,jdim
do i=1,idim
pt(i,j,l)=tmp(i,j,l)*factor
end doend doend doif( idout(1) /= 2008 .or. idout(2) /= 1 ) stop 11
ntime=(24*(idout(3)-idini(3))+(idout(4)-idini(4)))/12
if(ntime < 0 .or. ntime >= tmax ) stop 22
lrec=4*idim*jdim
kappa=rgas/cp

close(infile)finame='/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt'
foname='pt.gdt'

write(ioutfl,rec=1) pt

open(ioutfl,file=foname,iostat=ios,action='write', && recl=lrec*ldim,access='direct',form='unformatted')if(ios /= 0) stop 44

close(ioutfl)

gfortranの場合で、ifortの場合は4バイトを1単位としているので idim*jdim とする。

end program ptcalc

GrADSファイルの作成（温位データ）(続き)
GrADSデータ/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt から指定した日付の温位を計算
しGrADSデータファイルを作成する。
Cプログラム#include <stdio.h>#include <math.h>#define#define#definekappa = rgas/cp;
if( ( fpin = fopen( finame, "r" ) ) == NULL ) return 33;

IDIM 144
JDIM
73
LDIM
17

l = ( ntime*( 4*LDIM + 12 + 1 ) + LDIM )*IDIM*JDIM*sizeof(float);fseek( fpin, l, SEEK_SET );intidini[4]={2008,1,1,0},
idout[4]={2008,1,5,0},tmax=20;
float
plev[LDIM]={1000.0,925.0,850.0,700.0,600.0,
500.0,400.0,300.0,250.0,200.0,150.0,100.0,
70.0,50.0,30.0,20.0,10.0};
float
cp=1004.0,rgas=287.0,p0=1000.0;
float
tmp[LDIM][JDIM][IDIM],pt[LDIM][JDIM][IDIM];
char*finame="/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt",*foname="pt.gdt";

fread( tmp[0][0], sizeof(float), IDIM*JDIM*LDIM, fpin );for( l = 0; l < LDIM; l++ ) {factor = powf( p0/plev[l], kappa );for( j = 0; j < JDIM; j++ ) {for( i = 0; i < IDIM; i++ )pt[l][j][i] = tmp[l][j][i]*factor;}}

int main(){int
i,j,l,ntime;
float
kappa,factor;
FILE*fpin,*fpout;if( idout[0] != 2008 || idout[1] != 1 ) return 11;
ntime = ( 24*( idout[2] - idini[2] )+ ( idout[3] - idini[3] ) )/12;if( ntime < 0 || ntime >= tmax ) return 22;

fclose( fpin );if( ( fpout = fopen( foname, "w" ) ) == NULL ) return 44;
fwrite( pt[0][0], sizeof(float), IDIM*JDIM*LDIM, fpout );fclose( fpout );return 0;}

温位の緯度・高度断面表示

描画スクリプトの表示/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt から新たに作成した温位データのGrADSファイル
を使って、2008年1月5日00Zの東経140度に沿った緯度・高度断面図を下に示す。

図では熱帯で100hPa、中・高緯度で300hPaを境として、上下で安定度が大きく異なる。

成層圏

安定度の高い上部は成層圏、安定度の低い下部は対流圏に

前線面
対流圏

属する。
同じ対流圏でも、中・高緯度に比べ熱帯では等値線の間隔が開いており、大気の安定度がより低くなっていることが分かる。
図中の北緯20度の地上付近から北緯45度の300hPa付近にかけて等値線の密集した斜面が見られるが、これは前線面に対応している。

静力学平衡のチェック
総観規模現象のスケール（1000kmスケール）で見た場合、大気中では静力学平衡が良い近似で成り立っている。具体的には、気圧をp 、大気の密度をr 、重力加速度の大きさをg 、高度をz ととした場合

∂p
∂z

=− g

という関係が成立していることを意味する。この関係と大気の状態方程式

p= R Tを組み合わせると、以下の測高公式が求まる。但し、Φ はジオポテンシャル、T は気温、R は乾燥大気の気体定数である。

∂
∂p

RT
=−

または

∂
∂ ln pp=− R T

測高公式を気圧p1 からp2 まで（ p1 < p2 ）積分すると、p1 とp2 の間の層厚（thickness）が以下のように求まp2る。

 p1−  p2=∫ R T d ln p=R 〈 T 〉 ln p2/ p1 ⋯ 1
p1

ここで、<T >はp1 とp2 の間の平均気温で、次のように定義される。
p2〈 T 〉≡[∫ T d ln p]/ [ln  p2/ p1] ⋯ 2
p1

静力学平衡のチェック（GrADSで確かめてみよう）
1. /home/kisoen2/grads/sample_data/jra1st10dJAN08.ctl をGrADSで開く
2. 適当な緯度・経度の範囲と、時刻を「set」コマンドで設定する。
3. (1)式の左辺にあたる層厚を画像表示する。ただし、層厚は500hPaと850hPaの間とする。
注意：GrADSデータの「gph」は等圧面高度の値であるので、ジオポテンシャルの層厚は重力加速度g(=9.81) を掛ける必要がある。コマンドとしては

d 9.81*(gph(lev=500)-gph(lev=850))を入力する。
4. 表示した画像は画像ファイルとして保存する。
5. (1)式の右辺側を計算し画像表示する。(2)式で定義された平均気温を計算するのは面倒なので、<T>の代わりに500hPaと850hPaの中間にある700hPa気温で代用する。気体定数R(=287)とすると、表示コ
マンドはd 287*tmp(lev=700)*log(850/500)を入力する。
6. 低緯度の大気は大量の水蒸気を含んでいるので、乾燥大気の状態方程式は正しくない。水蒸気を含んだ湿潤大気では気温T の代わりに仮温度Tv を使えばよい。仮温度は比湿をq とすると

T v ≡T 10.608qと定義されているので、改めてコマンド

d 287*tmp(lev=700)*(1+0.608*sph(lev=700))*log(850/500)を入力して画像表示し、画像ファイルとして保存する。
7. 両者の図を比較してみる。

注意：図を比較するためには、両者とも「set cint」コマンドで等値線の間隔を揃えること。

静力学平衡のチェック（GrADSで確かめてみよう）（実例）
再解析データのGrADSデータ/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdtの2008年01月01日

00Zにおける等圧面高度と気温から静力学平衡のチェックをした例を下に示す。左側の図は500hPa高度と850hPa高度の差（層厚）に重力加速度を掛けてジオポテンシャルに変換した値（g*(z500-z850)）の分布で、単位はm2s-2である。一方、右側の図は500hPaと850hPaの間の平均気温の代わりに700hPaの気温を測高
公式に適用して推定した500hPaと850hPaのジオポテンシャルの差（R*Tv700*log(850/500)、Rは乾燥空気の気体定数、Tv700は700hPaの仮温度を表す）の分布である。層平均気温の変わりに700hPaの気温を使ったために両者の図には若干の違いはあるものの、静力学平衡が良く成り立っていることがわかる。

500hPaと850hPa高度から求めた層厚700hPa気温と測高公式から推定した層厚地衡風近似のチェック

差分近似の説明へ

大気の運動方程式は、水平方向の風速ベクトルをV 、コリオリパラメータをf 、ジオポテンシャルをΦ 、鉛直方向の単位ベクトルをk とすると、等圧面では

dV
dt=− f k ×V −∇ と記述されるが、中・高緯度では右辺の2項の絶対値が左辺のそれに比べ非常に大きいため、近似的に左辺を0とみなせる。この場合、風速ベクトルはジオポテンシャルで決まるが、この風を地衡風と呼ぶ。地衡風は中・高
緯度では実際の風の良い近似になっている。
東向きにx 座標、北向きにy 座標をとると、地衡風の東西成分ug と南北成分vg はそれぞれ

ug =−

1 ∂

f ∂y

, v g=

1 ∂

f ∂xと記述されるが、球面座標では地球の半径をa (=6.37×106m) 、経度をλ、緯度をφとすると

ug =−

1 ∂

f a ∂

, v g=1∂

f a cos ∂ と表せる。解析データの場合、それぞれの物理量が離散的な格子点上の値になっているので、微分は差分で次のように近似できる。

ug  i , j  ≃ −i , j 1− i , j −1
2 f j a  , v g  i ,  j  ≃i 1, j −i −1, j

2 f j a cos j  ⋯ 3

ここで、i , j はそれぞれ経度方向、緯度方向の格子点の番号で、東や北に向かって番号が増えるとする。また、
Dl, Dfはそれぞれ経度・緯度方向の格子間隔を表す。

地衡風近似のチェック（GrADSで確かめてみよう）
1. GrADSデータファイル/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt をFORTRANやC等のプログラミング言語で開いて、適当な時刻の500hPaにおける高度と風ベクトル（u, v）を読む。
2. 適当な緯度と経度の範囲（北極・南極と赤道は無意味なので避ける）を決めて、その範囲で(3)式に基づき地衡風を計算する。
3. 選択した緯度・経度の範囲の500hPa高度、元の風ベクトル、計算して求めた地衡風ベクトルを新たなGrADSファイルとして保存する。
4. 作成したGrADSファイルから、元の風ベクトルと地衡風ベクトルを比較する画像（500hPa高度も併せて表示させる）を作成する。
➢ (3)式の格子点のインデックスi, j に注意。元のGrADSデータファイ
ルは北から南に格子点値が配置されている。

備考
北➢ (3)式の格子点の間隔は、GrADSデータファイルでは
Dl=Df=2.5π/180ラジアンである。

Fi,j+1➢ (3)式のコリオリパラメータ fj は

Fi-1,jFi,jFi,j-1

Fi+1,jf j =2  sin j 

Df

で求められる。ここで、a (=6.37×106m) は地球の半径、
Ω (=7.292×10-5s-1) は地球の自転角速度を表す。

Dl
東

格子点値の配置➢ 風ベクトルを比較するときはベクトルの長さスケールを揃えること。
ベクトルの長さスケールを設定するGrADSコマンドは
set arrscl size magnitude
である（GrADSマニュアルを参照）。

地衡風近似のチェック（GrADSで確かめてみよう）（実例）
再解析データのGrADSデータ/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdtの2008年01月01日00Zにおける等圧面高度から地衡風を計算し実際の風と比較した例を下に示す。中緯度帯では風と地衡風は概ね良く一致しているが、30N,130Eから50N,180Eにかけてのジェット軸の入口・出口で大きなズレが見られる。この差はdV
dt=− f k ×V −∇ =−f k × V −V g という運動方程式から理解できる（ここで Vg は地衡風ベクトルを表す）。
GrADSファイルを作成したFORTRANのソース

500hPa高度と風ベクトル作成されたGrADSのコントロールファイル

GrADSファイルを作成したCのソース

500hPa高度と地衡風ベクトル

500hPa風と地衡風のベクトル差温度風のチェック
総観規模スケールで見ると、大気中では静力学平衡が非常に良い近似で成り立っている。

∂
∂p

=

−R T

p

∂
または

∂ ln pF ：ジオポテンシャル、T ：気温、R ：乾燥大気の気体定数=− R T一方、地衡風ug, vg はジオポテンシャルF とコリオリパラメータf から以下のように定義される。

ug=−

1∂

f ∂y

, vg =

1∂f ∂x

地衡風を気圧の対数で微分し、静力学平衡の式を代入すると、以下の式が得られる。

∂ ug∂ ln p

=

R ∂T
f ∂y

,

∂ vg
∂ ln p

=−

R ∂T
f ∂x

これらの式を気圧p1 からp2 （ p1 < p2 ）で積分すると、以下の温度風方程式が得られる。

ug  p1− ug  p2=−
ここで、<T> はR ∂〈 T 〉
f

∂y

ln

p2
p1

 , v g  p1− vg  p2=R ∂〈T 〉f∂x

p2

〈 T 〉≡[∫ T d ln p] / [ln  p2 / p1 ] ⋯ 5
p1で定義される、気圧p1 からp2 の間の平均気温を表す。

ln

p2
p1

 ⋯ 4温度風のチェック（続き）
温度風の定義：２つの等圧面間の地衡風のベクトル差（上層−下層）
温度風の意味（北半球の場合）

yy
寒v g  p1− vg  p20

寒Vg(p1)∂ 〈T 〉∂x0

Vg(p2)

暖

x

温
度
風 等
温
Vg(p2) 線

層平均気温が真東に向かって上昇しているとき、温度風は真北を向く。

温度風xVg(p1)∂ 〈T 〉∂y暖0ug  p1− ug  p20等温線

層平均気温が真北に向かって下降しているとき、温度風は真東を向く。

何れの場合も、北半球では温度風は層平均気温の暖気側を向かって右に見ながら、等温線に平行に吹く。

温度風のチェック（GrADSで確かめてみよう）
1. GrADSデータファイル/home/kisoen2/grads/sample_data/jra1st10dJAN08.gdt をFORTRANやC等のプログラミング言語で開いて、適当な時刻の500hPaと850hPa風ベクトル（u, v）と700hPaの気温T を読み込む。
2. 適当な緯度と経度の範囲（北極・南極と赤道は無意味なので避ける）を決めて、その範囲で(4)式に基づき温度風を計算する。但し、500hPaと850hPaの層平均気温<T > の代わりに両者の中間にある700hPa気
温T を使って温度風を計算する。その際、気温の偏微分は中央差分に置き換えて計算する（下記参照）。
3. 選択した緯度・経度の範囲の500hPaと850hPa風ベクトル（u, v）と700hPaの気温T、計算して求めた温度風ベクトルを新たなGrADSファイルとして保存する。
4. 作成したGrADSファイルから、元の風ベクトルから求めた500hPaと850hPaの風のベクトル差と、温度風ベクトルを比較する画像（温度風の計算に使用した700hPaの気温も併せて表示させる）を作成する。
5. 元のGrADSファイルで選択した緯度・経度の範囲の500hPaと850hPaの等圧面高度および両者の間の層厚の画像も作成し比較する。

北∂T

Ti,j+1

差分近似
Ti-1,jTi,jTi,j-1

Ti+1,j∂x≃

T i 1, j −T i −1, j
2a cos j   

,

∂T
∂y

≃

T i , j 1− T i , j −1
2a   6

a =6.37×10 m ,  = =2.5/180rad.

Df

Dl
東

注意：ジオポテンシャルは高度に重力加速度 g を掛けること。
JRAのGrADSデータは左図と異なり、北から南に並んでいる。

温度風のチェック（GrADSで確かめてみよう）（実例）
再解析データのGrADSデータ/home/kissoen2/grads/sample_data/jra1st10dJAN08.gdtの2008年01月01日00Z

における700hPa気温から温度風の式(4)に基づいて500hPa-850hPa間の温度風（地衡風のベクトル差、静力学
平衡と地衡風を仮定）を計算し、実際の500hPa-850hPa間の風のベクトル差と比較したもの。両者は概ね一致しているが、(4)式の500hPa-850hPa層間平均気温の代わりに700hPa気温で代用したことと、低緯度では地衡風近似が良くないために、25N以南の低緯度地帯とチベット高原のような高地（そもそも地表面が850hPaより上にある）では大きな不一致が見られる。
作成されたGrADSのコントロールファイル
GrADSファイルを作成したFORTRANのソース

700hPa気温と500hPa-850hPa
間の風のベクトル差

GrADSファイルを作成したCのソース

700hPa気温とそれから推定した
500hPa-850hPa間の温度風

500hPa高度（白線）、
850hPa高度（黒線）と両者の間の層厚（彩色）渦度の概念と性質●渦度の定義速度ベクトルをV としたとき、渦度ベクトルw はy

左図のように原点Oを中心として角速度Wで回転している
流れの場を考えると、中心からRの距離にある点での速度はRWで、方向は半径に直角だから、速度ベクトルのx、y 成分をそれぞれu、v とすると

RW
q

v

u=−R  sin =− yv =R  cos= xuWOだから、渦度ベクトルw のz 成分wz はxq≡∇ ×V で定義される。

R

z ≡

∂v
∂x

−

∂u
∂y

=2

になる。また、速度ベクトルのz 成分をw とすると w=0 だから

x ≡

∂w
∂y

−

∂v
∂z

=0 , y ≡

∂u
∂z

−

∂w
∂x

=0

よって、渦度は微小な領域での回転角速度の2倍の大きさで、回転軸の方向を向く、ことがわかる。

渦度の概念と性質（その２）●絶対渦度と相対渦度

W
w
v

u

半径がa で、角速度W で自転している地球表面を考える。ここでは球座標が定義できて、速度ベクトルV の経度l 方向の成分をu 、緯度f 方向の成分を
v 、半径（高さ）r 方向の成分をw とする。このとき渦度ベクトルw の高さ方向の成分wr は渦度の定義より

r ≡∇ ×V r =

a

1

[

∂v

a cos ∂ 

−

∂
∂

u cos]5という式で与えられる。
慣性系から見たときの速度Va=(ua, va, wa)は地球に相対的な（地上で観測した）速度V=(u, v, w)と、地球自体の回転速度Ve=(ue, ve, we)の和であるので、慣性系から見たとき、地上で観測したとき、地球自体の回転による、それぞれの渦度の高さ方向の成分をh、z、werとすると

fl

V a =V V e   ∇ ×V a r = ∇×V r  ∇×V e rから

= er が得られる。

地球自体の自転速度は V e = a  cos , 0,0 なので、（５）式を使うと

1
 er =

a cos

よって =  f

[

−

∂
∂

2

]

a  cos  =2 sin= f CoriolisParameter 

ここで、h は絶対渦度、z は相対渦度、f はコリオリパラメータである。

渦度の概念と性質（その３）●水平非発散場における渦度
東向きにx 座標、北向きにy 座標をとり、それぞれの風ベクトルの成分をu、v とする。鉛直p速度をw 、コリオリパラメータをf 、ジオポテンシャルをF とすると、大気の運動方程式は

∂u
∂t
∂v
∂t

u
u

∂u
∂x
∂v
∂x

v
v

∂u
∂y
∂v
∂y





∂u
∂p
∂v
∂p

=f v−

∂

⋯ 6

∂x

=− f u−

∂

⋯ 7

∂y

で表現される。(7)をx で偏微分したものから(6)をy で偏微分したものを差し引くと以下の渦度方程式が得られる。



∂
∂t

u

∂
∂x

v

∂
∂y



 f =− f 
∂v

∂u



∂u
∂x



∂v
∂y



−

∂
∂p





∂ u ∂∂p ∂y−∂v ∂ 
∂p ∂p



8

ここでz は相対渦度で ≡で定義される。また、コリオリパラメータf はx とt に依存しない−∂x ∂y
ことも利用した。
次に、総観規模現象における(8)式の各項の大きさを評価する。
総観規模現象では、代表的な水平スケールをL 、気圧のスケールをP 、水平速度のスケールをU 、鉛直P
速度のスケールをW とすると6と見積もれる。

5−1

L ≈10 m , P ≈10 Pa , W ≈10 Pa / s ⋯9

渦度の概念と性質（その４）(9)式のスケールを基に総観規模現象での(8)式の各項の大きさを見積もる。先ず、中緯度でのコリオリパ
−5 −1
−4 −1
ラメータの大きさは f ≈1×10 s
相対渦度の大きさは≈U / L ≈1×10 s
2であるから、(8)式左辺の大きさは

∂
∂t

,u

∂
∂x

U

≈

L

−10 −2s≈1×102程度になる。総観規模現象では発散は渦度の約1/10程度なので

∂u
∂x



∂v
∂y

≈0.1

U
L

−6 −1

≈1×10 sから右辺第1項は f と見積もれる。右辺第2項、第3項はそれぞれ



∂
∂p

≈

WU
PL

∂u ∂

−11 −2≈1×10s ,∂ p ∂y

,

∂v ∂

≈

∂ p ∂p



∂u
∂x

WU
PL



∂v
∂y



−10 −2≈1×10−11 −2≈1×10sというように左辺や右辺第1項に比べ1/10程度しかない。そこで渦度方程式は近似的に

もしくは




∂
∂t
∂
∂t

u

u

∂
∂x
∂
∂x

v

∂
∂y

v

と書ける。ここでh は絶対渦度を表す。

∂
∂y




 f ≃− f ≃− f 
=  f



∂u
∂x



∂u
∂x





∂v
∂y



∂v
∂y



⋯ 10

⋯ 11

s

渦度の概念と性質（その５）(11)式から、水平非発散の流れの場、つまり



∂
∂t

u

∂
∂x

v

∂
∂y



∂u
∂x



∂v
∂y

=0

が成り立つところでは

=0 ⋯ 12が成立する。これは、絶対渦度が水平の流れの中でラグランジュ的に保存されていることを意味する。

対流圏界面
発散

収束

∂u

500hPa

∂x



∂v
∂y

水平非発散

≃0

収束

発散
地

面

では、どの様な所で水平非発散になっているだろうか。
左図に示したように、対流圏全域にわたる
上昇（下降）流が発生したときは、対流圏界面付近で発散（収束）、地表面付近で収束（発散）という場になり、対流圏の中間にある500hPa付近では水平収束は０に近い。
従って、

500hPa付近では絶対渦度hが風の流れに沿って近似的に保存されている。

絶対渦度の保存を確かめよう●解析データから渦度を求める

絶対渦度h 、相対渦度z 、コリオリパラメータf の関係は

=  f ⋯ 13ui, j+1

相対渦度z は球座標系で表現すると

1
=

hi,j
vi-1, jDf
Dli , j =vi+1, jzi,jui, j-1

1

a cos j[vi 1, j −vi −1, j
2 

[

∂v

a cos ∂ 

−

∂

]

u cos  ⋯ 14∂となる。ここでa 、l 、f はそれぞれ地球の半径、経度、緯度をあらわす。
解析データは、緯度Df 、経度Dl の等間隔格子点上で定義されているので、相対渦度の格子点値は(14)式を差分化して−

 ui , j 1 cos j 1− ui , j −1 cos j −1
2 ]⋯ 15

から計算することが出来る。(15)式で添え字i は経度方向の格子点番号、添え字j は経度方向の格子
点番号を表し、i は西から東へ、j は南から北に向かって値が増えるとする。絶対渦度の格子点値は(13)式から以下の式で求められる。

 i , j = i , j  f j , f j =2 sin j ⋯ 16

絶対渦度の保存を確かめよう（その２）●作業の手順1. GrADSデータファイル/home/lessons/grads/sample_data/zuv500.gdtをFORTRANやC等のプログラミン
グ言語で開いて（データの構造は同じディレクトリのzuv500.ctlを参照）、全ての時刻（00Z23JAN2008～
18Z26JAN2008、6時間間隔、16個）の500hPaにおける風ベクトル（u, v）を読む。
2. 差分方程式(15)式と(16)式に基づき、全ての時刻の500hPaにおける絶対渦度h を計算し、そのGrADSファイル（データファイルとコントロールファイル）を作成する。ただし、GrADSデータファイルzuv500.gdtの格子点値の並びは、(15)式の添え字j の並びとは逆に北から南に並んでいることに注意すること。一方、東西方向の並びは(15)式の添え字i の並びと同じく西から東に並んでいる。(15)式では北極と南極の相対渦度が計算できないので、それぞれ北緯87.5度、南緯87.5度の相対渦度の平均値を極の渦度とする。
6

計算に際し、地球の半径 a =6.37×10 m 、格子点の間隔  = = / 72 を使うこと。
3. 作成した絶対渦度のGrADSファイルを使って、６時間毎の絶対渦度の分布図を画像ファイルにして保存し、WEBに図を表示させる。必要ならば、絶対渦度の図にzuv500.gdtを使って、ジオポテンシャル高度や風ベク
トルを付け加えてみてもよい。複数のGrADSを使って描画する方法はこちらを参照。
複数の画像をWEBで次々と表示させるには、HTML文書にJavaScriptのコントロールを追加すると便利である。その一例を kisoenshu2/abs_vor/ に示したので参考にすること。
4. WEB上に表示された絶対渦度の変化を眺め、絶対渦度の保存がどの程度当てはまっているのか確かめる。

絶対渦度の保存を確かめよう（その３）●JavaScriptを使ったHTML文書の例<html><head><meta http-equiv="content-type" content="text/html;
charset=utf-8"><title>絶対渦度の時間変化</title><style type="text/css"><!-body { margin:30px; background-color:navy; color:gold }--></style><script language="JavaScript"><!-var no,maxno=16;function showImg() {imgname="images/eta"+no+".png";
document.map.src=imgname;}function init() {no=1;
showImg();}function forward() {JavaScriptを
if( (++no) > maxno ) no=maxno;
使った制御部分
showImg();
ボタンクリックで}呼び出されるfunction backward() {if( (--no) < 1 ) no=1;
showImg();}function last() {no=maxno;
showImg();}--></script>

WEBページを開いたときに最初に呼ばれる関数<body onLoad="init()"><div align="center"><h1>絶対渦度の時間変化</h1><p><image src="images/eta1.png" name="map"><p><form name="control"><input type="button" value="|<" onClick="init()"><input type="button" value="<" onClick="backward()"><input type="button" value=">" onClick="forward()"><input type="button" value=">|" onClick="last()"></form><p>
備考：計算に使用したデータはJRA/JCDASの毎6時間解析値。等値線の間隔は1×10<sup>-5</sup>sec<sup>-1</sup></div></body></html>WEBページをコントロールするボタンの定義。
それぞれボタンが押されると、JavaScriptの関数、
init(), backward(), forward(), last()が呼ばれる。
関数は左側の「HEAD」部内で定義されている。
JavaScriptの変数「no」の値（1～16）に応じて画像ファイル「images/eta1.png」～「images/eta16.png」を表示させる。

JavaScriptを使うと、動的なWEBコンテンツを作ることができる。JavaScriptの詳細はインターネットで照会できる。

絶対渦度の保存を確かめよう（その４）●数値的なチェック
時刻t で座標r = ( x, y, p ) での水平風速をV = ( u, v )、絶対渦度を h(x,y,p,t) とすると[

  x ±u  t , y ±v  t , p, t ± t ≃  x , y , p, t =  x , y , p, t 



∂
∂t

u

∂
∂x

v

∂
∂y

∂
∂x

u

∂
∂y

v

∂
∂t

]

t

=0

が成り立つ。つまり、絶対渦度の保存が成り立っていれば、或る点の絶対渦度はその上流点の過去の値や、下流点の未来の値と一致する（ラグランジュ的に保存）。

●上流点・下流点の探査
絶対渦度がラグランジュ的に保存されているかど
うかチェックするためには上流点・下流点の探査が必要となる。或る緯度f 経度l の地点で風速ベクトルがV であるとき、Dt 前の時刻の上流点およびDt 後の時刻の下流点は、空間的にそれぞれベクトルV Dt 進んだ、もしくは戻った点を探せば良い。
但し、球面上であるため、探査は風速ベクトルVに接する大円上で行う必要がある。

下流
V 点

f

上流
点

l

大円

絶対渦度の保存を確かめよう（その５）●探査すべき大円の決定
緯度f 、経度l を通る大円は(17),cos cos =cos cos −sin  cos sin ⋯ 17

(18), (19)式で定義出来る。ここで、sin  cos=sin  cos cos  cos sin ⋯ 18

a 、b は決めるべきパラメータ、g は

sin=sin  sin 

大円上で-pから+pの値をとる。
この大円は緯度f 、経度l で⋯ 19

 = 0 となるとする。この点で、風速ベクトルV が緯線と角度d を持つとき、風速ベクトルに接する大円は以下の式から唯一に決定出来る（つまりパラメータa 、b が決まる）。

cos =cos  cos ⋯ 20 , sin  0=

sin 
sin 

⋯ 21 , cos− =cos 0cos

なお、 =± / 2 のときは

= , = / 2 ,  0=

if = / 2 ⋯ 23

=− ,  =/ 2 ,  0=−
から a 、 b 、 g0 が決定される。

if =− / 2 ⋯ 24⋯ 22

絶対渦度の保存を確かめよう（その６）地球の半径をa とすると、上流点・下流点のg の値g- 、g+ は

 -= 0−∣V ∣ t ,  = ∣V ∣ t ⋯ 23
+

a

0

a

から求められ、それらに対応した緯度・経度が(17), (18), (19)式から得られる。

●上流点・下流点での空間内挿求められた上流点・下流点は必ずしも格子点上にないため空間内挿が必要となる。右図のように上流点の緯度・経度をそれぞれf- 、l- とすると、上流点

大円

hi, j+1

下の式で内挿して求める。

f-fhi+1, j+1 j+1
上流点  r − V  t , t − t hi, jhi+1, jl i l-

li+1の絶対渦度は上流点を囲む４つの格子点上の t-Dt における絶対渦度から

 i 1−-  j 1−- t − t= i , jt − t i , j 1

i 1−i  j 1− j 
 i 1−-  -− j 
i 1−i  j 1− j 

t − t i 1, j

 -−i  j 1− - 
i 1−i  j 1− j 

t − t i 1, j 1

- −i - − j 
i 1− i  j 1− j fj絶対渦度の保存を確かめよう（その７）00Z25AUG2008の解析値計算に使用したデータはJRA/JCDASの500hPa風ベクトル
00Z23JAN2008〜18Z26JAN2008の値/home/lessons/grads/
sample_data/zuv500.gdt

12時間前の上流点から推定した値

12時間後の下流点から推定した値

可視化ツールを使った解析データの利用と、気象を支配する物理法則の体験的理解

最終スライド