# draw 1OOO-day averages of the zonal mean potential temperature
# (for dry atmosphere and flat topography)
#
# PotTemp = Temp * sigma ** (- R/C_p)
# 
# 更新履歴
# 
# 2009/03/29  納多 哲史   新規作成
# 

require "numru/ggraph"
include NumRu

# 定数の設定
# 乾燥空気の気体定数
r = 2.87e+2     # [J K-1 kg-1]
# 乾燥空気の定圧比熱
cp = 1.004e+3   # [J K-1 kg-1]

# 引数のチェック
if ARGV[0] == "-h" or ARGV[0] == "--help" then
  help()
  exit 0
elsif ARGV.length != 2 then
  puts "ERROR: The number of argument is invalid."
  help()
  exit 1
end

file = ARGV[0]
varname = ARGV[1]

if !FileTest.file?(file) then
  puts "ERROR: No such file."
  exit 2
end

gphys = GPhys::IO.open(file, varname)
gphys = gphys.cut('time'=>200..1200)

# 以下の方法はかなり恰好悪い. もっと簡単に重み付けをする方法があるはず.

nlon = gphys.coord('lon').length
nlat = gphys.coord('lat').length
nsig = gphys.coord('sig').length

weight = NArray.sfloat(nlon, nlat, nsig)
wsig = gphys.coord('sig').to_a

# 温度に掛ける係数を計算. 地面が平坦であることを仮定.
for k in 0...(nsig-1)
  weight[0..(nlon-1), 0..(nlat-1), k] = wsig[k] ** (-r/cp)
end

# 温度から温位に変換
gphys = gphys * weight

gphys = gphys.mean('lon', 'time')

DCL.gropn(2)

GGraph.set_fig('viewport'=>[0.15,0.80,0.15,0.6])

DCL.sgpset('lfull',true)     # 全画面表示
DCL.sgpset('lcntl', false) ; DCL.uzfact(0.6)

GGraph.set_linear_contour_options('interval'=>5, 'max'=>500)

# 最上層より上はコンターが混みすぎて, DCL でエラーになるため表示しない
GGraph.contour(gphys.cut('sig'=>0.025..1))

DCL.grcls