（クリックでフルサイズ表示）
Ruby と GPhys を使いこなせば，地球流体の研究で行う様々なデータ解析ができます． 本章はそのための入門編です．この先，GPhysレシピ を参考にさらに経験を積めばなんでも (^_^) できるようになります． *1
目次
ある座標軸に沿って平均するという操作は，これまでも出てきました． ここではより包括的に解説します． mean は，格子点に関するアンサンブル平均です （格子間隔にもとづく重み等はかけません）． その仕様は次のようになってます．
GPhys#mean(*dims)
ここでは Ruby 流の記述をしました．これにも慣れましょう：
mean と同じ仕様のメソッドには次があります：
いずれも NArray, NArrayMiss のメソッドです．データに欠損があるときは NArrayMiss が使われます． NArrayMiss は，欠損は飛ばして処理します（平均するときの母数にもカウントしません）． GPhys は処理の中身はこれら委ねます．
では，サンプルプログラムです．
ソースコード sample_mean_etc.rb:
require "numru/gphys"
include NumRu
gp = GPhys::IO.open "air.2012-01.nc","air"
gp = gp.cut("level"=>925..20, "lat"=>40)[true,{0..-1=>2},true]
#^サブセット指定 (この時点ではまだ読み込まれない)
print "gp :\n"
p gp
print "pressure levels:\n"
p gp.coord("level").val
print "\ngp.mean(0) :\n"
p gp.mean(0)
print "\ngp.mean('lon','time') :\n"
p gp.mean('lon','time')
print "\ngp.mean :\n"
mn = gp.mean
p mn.class, mn
print "\nmean,max,min,stddev values:\n"
p gp.mean('lon','time').val, gp.max('lon','time').val,
gp.min('lon','time').val, gp.stddev('lon','time').val実行してみましょう．
$ ruby sample_mean_etc.rb
gp :
<NumRu::GPhys grid=<3D grid <axis pos=<'lon' shape=[144] subset of a NumRu::VArrayNetCDF>>
<axis pos=<'level' shape=[8] subset of a NumRu::VArrayNetCDF>>
<axis pos=<'time' shape=[31] subset of a NumRu::VArrayNetCDF>>>
data=<'air' shape=[144, 8, 31] subset of a NumRu::VArrayNetCDF>>
pressure levels:
NArray.sfloat(8):
[ 925.0, 700.0, 500.0, 300.0, 200.0, 100.0, 50.0, 20.0 ]
gp.mean(0) :
<NumRu::GPhys grid=<2D grid <axis pos=<'level' shape=[8] subset of a NumRu::VArrayNetCDF>>
<axis pos=<'time' shape=[31] subset of a NumRu::VArrayNetCDF>>>
data=<'air' sfloat[8, 31] val=[276.323852539062,266.1376953125,250.100402832031,226.781204223633,...]>>
gp.mean('lon','time') :
<NumRu::GPhys grid=<1D grid <axis pos=<'level' shape=[8] subset of a NumRu::VArrayNetCDF>>>
data=<'air' sfloat[8] val=[274.950927734375,264.612152099609,248.710037231445,224.867828369141,...]>>
gp.mean :
NumRu::UNumeric
234.546216957885 degK
mean,max,min,stddev values:
NArrayMiss.sfloat(8):
[ 274.951, 264.612, 248.71, 224.868, 217.77, 214.26, 213.517, 217.683 ]
NArrayMiss.sfloat(8):
[ 287.55, 278.02, 261.12, 235.8, 236.57, 226.35, 226.35, 235.0 ]
NArrayMiss.sfloat(8):
[ 253.87, 245.42, 230.77, 215.32, 205.35, 202.62, 200.57, 200.32 ]
NArrayMiss.sfloat(8):
[ 6.73178, 6.87236, 6.20152, 4.12792, 6.29422, 4.00543, 4.48009, 6.30728 ]
ソースの最初は require "numru/gphys" となってます．可視化を行わない場合，numru/ggraph でなく numru/gphys を require することで起動が少し早くなります*2．
プログラムでは，最初に gp.cut.... の行でサブセットをとります．
GPhys のサブセット切り出し（[]やcutによる）は，
もとのデータのサブセットへのマッピング法を格納したオブジェクトを作るだけなので，
この時点ではデータは読み込みません（だからもしもサブセットに書き込めば参照元に書き込まれます）．
このことは，出力の最初のほうに
subset of a NumRu::VArrayNetCDF と現れることから伺えます．
データ読み込みは gp.mean(0) などの演算を行う際にはじめて行われます．
その際は必要なサブセットだけが読み込みまれます．
ですから，サブセット切り出しは演算の前に行うのがよいです．
結果をみると，mean の戻り値は引数ありの場合 GPhys, なしの場合 UNumeric であることがわかります（省略せずに書くと NumRu::GPhys と NumRu::UNumeric です）．gp.mean(0) の結果の表示では， data=<'air' sfloat[7] val=[274.950927734375,264.612152099609,... となってることに注目してください．演算の結果はメモリー上に NArray または NArrayMiss として置かれます．つまり，同じ GPhys オブジェクトでも， 最初はデータはファイル中にありますが，演算の結果できたものは全体が計算機のメモリー上に保持されるのです． このように違っていても同じ GPhys なので同じように描画できるというのは， これまでの章で体験しました（こういうことが簡単に実現できるのがオブジェクト指向言語の利点です）． なお，演算結果をファイルに書き出す方法は後で示します．
プログラムの最後では，緯度40度の各気圧面から得られたいくつかの統計量を表示します． 母集団は経度および時間の2次元断面にかかる 144×31=4464 個のデータです． 入力データは欠損があり得る形で NetCDF ファイルに格納されているので 結果は NArrayMiss クラスですが，実際には欠損はありません （はじめる前にの末尾を参照）．
次のプログラムは2012年1月後半の500hPa等圧面高度の半旬（5日）平均値（コンター）と， その気候値からの偏差（色）を表示します． 最初の図（1月16-20日）にみられる日付変更線付近の高緯度の高気圧偏差はブロッキング高気圧です． タイトルに平均期間が表示されるよう工夫してあることに注意してください．
ソースコード sample_pentad.rb:
# coding: cp932
require "numru/ggraph"
include NumRu
plev = 500
z = GPhys::IO.open('hgt.2012-01.nc', 'hgt').cut('level'=>plev)
zc = GPhys::IO.open('ncep2.Jan.clim.1981-2010.nc', 'hgt').cut('level'=>plev)
iws= 1
DCL.sgscmn(4) # カラーマップ番号
DCL.gropn(iws) # iws : DCL出力デバイス．1,4:画面, 2:PS
DCL.sgpset('isub', 96) # 下付き添字制御文字を '_' から '`' に
DCL.glpset('lmiss',true) # DCLの欠損値処理を on に
GGraph.set_fig "itr"=>32,"viewport"=>[0.15,0.85,0.15,0.85] # 32:正距方位図法
GGraph.set_map "coast_world"=>true
(16..26).step(5) do |d0| # 2012年1月の日付でループ
t0 = DateTime.new(2012,1,d0) # 2012-01-d0 00UTC (時刻指定省略のため)
t1 = t0 + 4.9 # 5日間の終わり(データは毎正時なので +4 でも可)
zp = z.cut("time"=>t0..t1).mean("time") # 半旬(5日)平均
zpd = zp - zc # 気候値からのずれ
GGraph.tone zpd,true,"max"=>500,"min"=>-500,
"annot"=>false, # 右脇のアノテーション非表示
"title"=>"Z#{plev} #{t0.strftime('%Y/%m/%d')}-#{t1.strftime('%d')}"
GGraph.contour zp,false,"int"=>100
GGraph.color_bar
end
DCL.grcls結果の一枚：
（クリックでフルサイズ表示）
データ処理した結果をファイルにとっておきたいことは多いでしょう． GPhys では次のように簡単に NetCDF ファイルに書き出せます．
ソースコード sample_output.rb:
require "numru/gphys" include NumRu gp = GPhys::IO.open "air.2012-01.nc","air" gpmn = gp.mean(-1) ofile = NetCDF.create "air.2012-01.monmean.nc" GPhys::IO.write ofile, gpmn ofile.close
実行と結果確認：
$ ruby sample_output.rb
[horinout@localhost trunk]$ ncdump -h air.2012-01.monmean.nc
netcdf air.2012-01.monmean {
dimensions:
lon = 144 ;
lat = 73 ;
level = 17 ;
variables:
float lon(lon) ;
...中略...
float air(level, lat, lon) ;
...後略...
出力のポイントは３つ：
NetCDF.create(パス)（新規）または，NetCDF.open(パス,"a")（追記）でファイルをオープンするGPhys::IO.write(ファイルオブジェクト，GPhysオブジェクト) で出力．ファイルオブジェクト.close で終了処理です．GPhys::IO.write は，座標軸を含む GPhys オブジェクトを構成する要素をすべてファイルに書き出すので，これだけで自己記述性のあるデータが書き出されます．
なお，ファイル中の名前は重複しないようにしなければなりません． よって，一つのファイルに同じ名前の複数の変数（例えば "air" の経度平均と時間平均） を出力したければ少なくともいずれかは GPhys"name= メソッドを使って改名する必要があります．
GPhys には他にもいろいろな数値計算手法， さらに気象学・地球流体力学的なデータ解析手法を提供するライブラリがあります． 数値微分，フーリエ解析，球面幾何，気象学的な諸量（温位，渦位等）の計算等々です． それらについては， GPhysのマニュアル のほか， GPhysレシピ でも紹介する予定です． 以下では．GPhys付属のライブラリ利用例を示し，自分でライブラリを作る例を示します．
次の例は，渦位（ポテンシャル渦度; potential vorticity; PV）を計算して図示します．
ソースコード sample_pv.rb:
# coding: cp932
# 実行方法：
# $ ruby sample_pv.rb [温位 [アニメーション間隔 [最初の日付 [最後の日付]]]]
# 実行例：
# $ ruby sample_pv.rb
# $ ruby sample_pv.rb 700
# $ ruby sample_pv.rb 320 2
# $ ruby sample_pv.rb 320 0 14 20
#
require "numru/ganalysis"
require "numru/ggraph"
include NumRu
#< コマンドラインオプション >
theta = ( ARGV[0] || 320.0 ).to_f # 温位面を指定
tsleep = [ ( ARGV[1] || 1 ).to_f, 0.0 ].max # 第2引数は描画間隔(秒)
fday = ( ARGV[2] || 10 ).to_i # 第3引数は最初の日付(1月何日か)
lday = ( ARGV[3] || 22 ).to_i # 第4引数は最後の日付(1月何日か)
wait = ( tsleep <= 0.0 ) #->true/false; 0以下(true)ならマウスクリックを待つ
theta_levs = [theta]
#< 決め打ちの定数 >
iws = 1
#< 使用データ >
cut = {"lat"=>0..90, # 表示する範囲にあわせて制限しておくと無駄な計算が減る
# "level"=>1000..100, # 高度制限（計算可能な温位の範囲は減るが）
"time"=>Date.new(2012,1,fday)..Date.new(2012,1,lday)}
tmp = GPhys::IO.open('air.2012-01.nc', 'air').cut(cut)
u = GPhys::IO.open('uwnd.2012-01.nc', 'uwnd').cut(cut)
v = GPhys::IO.open('vwnd.2012-01.nc', 'vwnd').cut(cut)
theta = GAnalysis::Met.temp2theta(tmp) # 温位をもとめる
pv = GAnalysis::Met.pv_on_theta(u, v, theta, theta_levs) # PVをもとめる
#< 描画準備 >
DCL.sgscmn(10) # カラーマップ番号指定
DCL.swpset("iwidth",500) # 画面の幅
DCL.swpset("iheight",500) # 画面の高さ
##DCL.swpset("ldump",true) # 画像ファイルダンプ
DCL.swpset("lwait",wait) # 次の描画の前にマウスクリックを待つ
DCL.swpset("lalt",true) # 裏で描画（パラパラアニメ用）
DCL.gropn(iws) # DCL描画装置初期化 (iws=1,2,or 4)
DCL.sgpset('isub', 96) # 下付き添字を表す制御文字を '_' から '`' に
# (アンダースコアはよく使われるので，エラー防止のため)
DCL.glpset('lmiss',true) # DCLの欠損値処理を on に．
GGraph.set_fig "itr"=>32,"viewport"=>[0.05,0.81,0.12,0.88] # 32:正距方位図法
GGraph.set_map "coast_world"=>true
#< 描画 >
nt = pv.shape[-1]
for it in 0...nt
GGraph.tone pv[false,it].cut("lat"=>0..87.5),true, # 極を除いて描画
"keep"=>(it!=0) # "keep"=>true で前回のカラーマップをそのまま使う
GGraph.color_bar
sleep tsleep
end
DCL.grcls最初の行は，これまでのソースにはない
require "numru/ganalysis" となっています．
こうすると GPhys において応用的なライブラリを集めた GAnalysis
というモジュールが利用できるようになります．
GAnalysis は GPhys パッケージの一部ですが，
このように陽に指定しない限り読み込まれないようになっています．
これにより，GAnalysis が必要ない場合の GPhys
の立ち上げが速くなってます．
上のプログラムでは渦位は
pv = GAnalysis::Met.pv_on_theta(u, v, theta, theta_levs)
によって計算されます．これは GAnalysis 中の
Met というモジュール中の pv_on_theta
というメソッドで，気圧座標の入力データを元に等温位面（一般に複数個）上での渦位を計算します
（GAnalysis::Metのオンラインマニュアル参照）．
４つの引数は順番に，東西風，南北風，温位（以上同じ格子点上で定義された
GPhys オブジェクト），渦位を求めたい温位面（数値のArray）です．
プログラム冒頭に書かれているように，４つのコマンドラインオプション（引数）があります． 最初の引数では温位面を指定します．以下はデフォルトの 320 K での結果例です． 上述の sample_pentad.rb で表示したブロッキング高気圧の形成過程を示します． 低緯度の低渦位の空気塊が高緯度に伸び切離したことがわかります．
実行結果より：
さて，付属のライブラリがどんなに充実してもすべての研究者の需要はみたせません． データ解析手法を自分で実装しライブラリ化して使いまわせるようになりましょう． ライブラリ化といっても Ruby では簡単で，独立したファイルにするだけです． そのファイルを require すれば利用できます． （require については，一般的な Ruby の参考書等を参照してください）． 以下に作成例をひとつ示します．サンプルプログラム付きのライブラリになります．
とりあげるのは Takaya&Nakamura (2001)の波活動度フラックスです． それは次の式で表されます (Takaya and Nakamura, 2001, J.Atmos.Sci, 608-627, Eq.(32))．
ここで，U = [U,V,0] は「基本場」の東西，南北流，ψ' は準地衡流線関数の偏差，p は気圧を 1000 hPa で規格化したものです． x,yは東西，南北で球面の効果を無視しています （球面版の実装も難しくありませんが，ここでは簡単のため x,y 版にします）． Cu は波の位相速度に関連する量で停滞波ではゼロになるので，ここでは Cu M の項は無視します（いつでも無視できるわけではないので適用限界があります）． なお，簡単のため，水平成分だけ計算することにします．
この波活動度フラックスの計算を実装してみます． （慣れればこの程度の複雑さのものはすぐ作れます．しかも GPhys メソッドとして作るため入力に関しかなり汎用になります．） では，ちょっと長いですがソースをざっと眺めてみましょう． ライブラリの主な（実質的に唯一の）部分はメソッド tn2001fluxCSH です． 中をみると，数式との対応が見て取れるのではないでしょうか． （なお，上記の事情を勘案して背景風の取り扱いについて二通りの計算が用意してあります．）
ソースコード tn2001flux.rb:
# coding: cp932
require "numru/gphys"
require "numru/gphys/derivative"
module NumRu
class GPhys
# Wave activity flux by Takata&Nakamura(2001) (Cartesian; Stationary; Horizontal component only)
#
# Simplified version neglecting the C_U*M term, which is zero for stationary
# waves. (if the waves are not stationary, the result is the wave activity
# flux relative to their phase propagation)
#
# HISTORY
# * T.Horinouchi 2012-12-07
#
# CAUTION: The above neglect is equivalent to assume that the
# ground-based phase speed is much smaller than the backgound wind
# speed, WHICH MAY NOT BE VALID IF THE BACKGROUND WIND IS WEAK!!
#
# self must be geopotential height anomaly (usu. deviation from
# climatology) on isobaric surface(s)
# with longitude and latitude as the 1st & 2nd coordinates
#
# ARGUMENTS:
# * ubk (GPhys/nil) : background zonal wind (w/ same grid as self).
# If nil (special case), bg wind is not used (equiv to assume vbk==0)
# * vbk (GPhys/nil) : background merid wind (w/ same grid as self).
# If nil (special case), bg wind is not used (equiv to assume vbk==0)
# * prs (nil or UNumeric; default: nil) : pressure level (good for single
# pressure-level data); if nil, the 2nd dim is assumed to be pressue
# * reflat (Numeric; default: 45degN) : referential latitude (in degN)
#
# RETURN VALUE:
# * [wx ,wy] : where wx&wy are the x&y components of the
# wave activity flux (the vertical component is not derived)
#
def tn2001fluxCSH(ubk, vbk, prs=nil, reflat=45.0)
#< assumptions on dimensions : you can modify it to configure it! >
xd = 0 # longitudinal dimension
yd = 1 # latitudinal dimension
pd = 2 # [used only if prs is nil] pressure dimension
#< constants >
g = UNumeric[9.8, "m.s-2"] # gravity
a = UNumeric[6.37e6, "m"] # Earth's radius
trot = 86400.0 # Earth's rotation in sec (to calc f0 below)
y0 = Math::PI * reflat / 180.0 # reference latitude (radian)
f0 = UNumeric[4*Math::PI / trot * Math::sin(y0), "s-1"] # Coriolis para
#< preparation >
za = self # geoptential height anomaly
psi = za * (g/f0) # quasi-geostrophic stream function
x = psi.axis(xd).pos.convert_units("rad") * (a*Math::cos(y0))
y = psi.axis(yd).pos.convert_units("rad") * a
x.units = "m" # rad.m -> m
y.units = "m" # rad.m -> m
psi.axis(xd).set_pos(x) # replace lon ax for differentiation
psi.axis(yd).set_pos(y) # replace lat ax for differentiation
prs = psi.axis(pd).to_gphys if !prs
#< calculation >
psi_x = psi.threepoint_O2nd_deriv(xd)
psi_xx = psi.deriv2nd(xd)
psi_y = psi.threepoint_O2nd_deriv(yd)
psi_yy = psi.deriv2nd(yd)
psi_xy = psi_x.threepoint_O2nd_deriv(yd)
if ubk && vbk
#Case: TN2001 Eq.32 except the CuM term
c = prs / (2*UNumeric[1000.0,"hPa"]) / (ubk**2 + vbk**2).sqrt #rho0/2U
ub = ubk*c
vb = vbk*c
wx = ub*(psi_x**2 - psi*psi_xx) + vb*(psi_x*psi_y - psi*psi_xy)
wy = ub*(psi_x*psi_y - psi*psi_xy) + vb*(psi_y**2 - psi*psi_yy)
wx.long_name = "WAF-x (stationary)"
wy.long_name = "WAF-y (stationary)"
else
#Case: Further simplification assuming V=0 (or the zonal momentum flux)
c = prs / (2*UNumeric[1000.0,"hPa"])
wx = (psi_x**2 - psi*psi_xx) * c
wy = (psi_x*psi_y - psi*psi_xy) * c
wx.long_name = "WAF-x (stationary,V=0)"
wy.long_name = "WAF-y (stationary,V=0)"
end
#< post processing >
grid = za.grid
wx = GPhys.new(grid,wx.data) # use the original grid
wy = GPhys.new(grid,wy.data) # use the original grid
wx.name = "wx"
wy.name = "wy"
[wx,wy]
end
end
end
######################################################################
## Main program (test part)
if __FILE__ == $0
require "numru/ggraph"
include NumRu
def prep_grph
iws= 1
DCL.sgscmn(4) # カラーマップ番号
DCL.swpset("ldump",true) if iws==4 # 画像ファイルダンプ
DCL.swpset('iwidth',1000)
DCL.swpset('iheight',500)
DCL.gropn(iws) # iws : DCL出力デバイス．1,4:画面, 2:PS
DCL.sldiv("y",2,1)
#DCL.sgpset("lfull",true)
DCL.sgpset('isub', 96) # 下付き添字制御文字を '_' から '`' に
DCL.glpset('lmiss',true) # DCLの欠損値処理を on に
#GGraph.set_fig "itr"=>10,"viewport"=>[0.05,0.85,0.2,0.5] # 10:経度緯度座標
GGraph.set_fig "itr"=>32,"viewport"=>[0.06,0.86,0.1,0.9] # 32:正距方位図法
GGraph.set_map "coast_world"=>true
end
plev = 300
latran = 20..75
z = GPhys::IO.open('hgt.2012-01.nc', 'hgt').cut('level'=>plev,"lat"=>latran)
clifl = 'ncep2.Jan.clim.1981-2010.nc'
zcli = GPhys::IO.open(clifl,'hgt').cut('level'=>plev,"lat"=>latran)
ubk = GPhys::IO.open(clifl,'uwnd').cut('level'=>plev,"lat"=>latran)
vbk = GPhys::IO.open(clifl,'vwnd').cut('level'=>plev,"lat"=>latran)
prep_grph
(8..18).step(5) do |d0| # 2012年1月の日付でループ
t0 = DateTime.new(2012,1,d0) # 2012-01-d0 00UTC (時刻指定省略のため)
t1 = t0 + 9.9 # 5日間の終わり(データは毎正時なので +4 でも可)
zp = z.cut("time"=>t0..t1).mean("time") # 半旬(5日)平均
za = zp - zcli
[ [ubk, vbk], [nil,nil ] ].each do |ub,vb|
wx,wy = za.tn2001fluxCSH(ub, vb, UNumeric[plev,"hPa"], 45.0 )
su = !ub.nil? ? "with UV" : "w/o UV"
GGraph.tone za,true,"max"=>500,"min"=>-500,
"annot"=>false, # 右脇のアノテーション非表示
"title"=>"Z#{plev} & WAF #{su} #{t0.strftime('%Y/%m/%d')}-#{t1.strftime('%d')}"
GGraph.contour zp,false,"int"=>100
GGraph.vector wx,wy,false,"xint"=>3,"yint"=>2,"unit"=>true,"fact"=>1.1
GGraph.color_bar "vcent"=>0.4,"vlength"=>0.3
end
end
GGraph.tone ubk,true,"max"=>75,"min"=>-75,"int"=>15,"title"=>"clim U & (U,V)"
GGraph.vector ubk,vbk,false,"xint"=>3,"yint"=>2,"unit"=>true,"fact"=>2
GGraph.color_bar "vcent"=>0.4,"vlength"=>0.3
uu = (ubk**2 + vbk**2).sqrt
GGraph.tone uu,true,"max"=>70,"min"=>0,"int"=>5,"title"=>"clim WindSpeed & (U,V)"
GGraph.vector ubk,vbk,false,"xint"=>3,"yint"=>2,"unit"=>true,"fact"=>2
GGraph.color_bar "vcent"=>0.4,"vlength"=>0.3
DCL.grcls
endメソッド tn2001fluxCSH は
module NumRu class GPhys end end
にくるまれてるので，GPhys のメソッドとして定義されていることがわかります． 微分は差分で近似しました．経度，緯度座標はあらかじめ x = a cosφ0 λ, y = aφ に変換しておきます（aは地球半径，λは経度，φは緯度，φ0は基準となる緯度です）． 使った差分法は２次精度で不等間隔格子に対応します（一階微分は等間隔格子なら中央差分になります）． 仕様についてはマニュアル(Reference manuals of the Derivative library) を参照してください．
ライブラリ定義の下の部分()は
if __FILE__ == $0 end
で囲われています．このファイルを他のファイルから require するとこの中は実行されず，このファイルを直接 ruby コマンドの引数とした場合のみ実行されます． ライブラリの利用法のサンプル兼テストプログラムになっている次第です． このやり方は，Ruby でライブラリを作る時の常道の一つです *3．
では実行してみましょう．
% ruby tn2001flux.rb
結果の一枚：
（クリックでフルサイズ表示）
参考：背景風（最終ページ）：
（クリックでフルサイズ表示）
ここでは背景場は気候値にしました． 本計算には停滞性の仮定があるため，先程より長い10日間の平均をとって， 気候値からのずれに対し波活動度フラックスを計算しました． 左右の図の違いは矢印で表した波活動度フラックスの計算法のみです． 左は Cu = 0 とした式(32)そのままですが， 右は背景南北風がゼロと仮定した場合に相当します． 左右が大きく異なるところは，|V| が |U| に比して小さいと言えないわけですが， 一般に |V| は小さいので，背景風自体が小さい可能性があります（上図のブロッキング高気圧あたりがそうです）． そうなると，背景風に比して位相伝播を無視するという停滞性の仮定はよくないかもしれないので注意が必要です（なお，今回無視した CuM を計算するには波の位相速度を特定しないとなりません）． ちなみに，今回のように CuM を無視すると，波の位相伝播から相対的にみた波活動度のフラックスになります．その意味で，停滞性の波でなくても意味がある計算ではあります．
なお，本例に対し論文の著者の高谷康太郎さんより次のようなコメントをいただきました．
”アラスカ付近のブロッキングを対象とされていますが、例としては一番解釈の難しいケースを取り上げられたな、という感じはします。御指摘の通り、アラスカ付近は気候値でも背景場が歪んでおり背景風速の弱い領域です。...(中略)...チュートリアルとしては、むしろ、ジェット中のブロッキングであるヨーロッパ付近やユーラシア大陸上のものの方が良かったかも知れません。もちろん、ブロッキング高気圧そのものは非線形現象ですが、ブロッキングへの上流からのロスビー波の入射、またはブロッキングから下流への射出、などは線形性の強い現象であり、フラックスの例としてはより良いかと思います(TN2001またはTakaya and Nakamura 2005など)。”
*1「なんでも」というのは言いすぎかもしれませんが，GrADSなどと違い Ruby
は本格的なプログラミング言語なので自由度は非常に高いです．
計算速度がネックになれば C で拡張することで補えますし．
*2numru/ggraph を
require するとその中で numru/gphys が require されます
*3もっと厳格に「テスト駆動開発」もできますが，ここでは一番簡単なやり方として紹介します．いずれにしてもサンプルは欲しいですし．
