# -*- coding: euc-jp -*-
#
# 2011/07/13 drawgm.rb を元に作成
#


require "numru/ggraph"
include NumRu


dir    = '/home2/yukai-home2/NCEP'
vname1 = 'sst_T62'
vname2 = 'ulwrf'
vname3 = 'uswrf'
vname4 = 'dswrf'

file2 = 'ulwrftop.mon.mean'
file3 = 'uswrftop.mon.mean'
file4 = 'dswrftop.mon.mean'


#
# NECP の データは 1988-2007 年度の各月平均を使用する
#
ts2 = 17417760
te2 = 17592240

#
# 緯度方向にかける重みの読み込み
#
gphysw = GPhys::NetCDF_IO.open(".."+'/'+vname1+".nc", 'lat_weight')
weight = gphysw.val
weight = weight / weight.sum(0)

#p weight.sum(0)
#p weight.reshape!(1,gphysw.shape[0])


#
# 各 nc ファイルの読み込み
#
gphys2 = GPhys::NetCDF_IO.open(dir+'/'+file2+".nc", vname2)
gphys3 = GPhys::NetCDF_IO.open(dir+'/'+file3+".nc", vname3)
gphys4 = GPhys::NetCDF_IO.open(dir+'/'+file4+".nc", vname4)

gphys5 = gphys4 - gphys3 # 射出太陽放射
gphys5.name      = 'outgoing shortwave'
gphys5.long_name = 'outgoing shortwave'

#
# 変数 time の定義
#
time = gphys2.coord('time').val


gphys2 = ( gphys2 * weight.reshape!(1,gphysw.shape[0]) ).sum(1).mean(0).cut('time'=>ts2..te2).mean('time')
gphys5 = ( gphys5 * weight.reshape!(1,gphysw.shape[0]) ).sum(1).mean(0).cut('time'=>ts2..te2).mean('time')

p 'OLR' + ': ' + gphys2.to_s
p 'OSR' + ': ' + gphys5.to_s

fileunit = File.open("radiation-budget.txt", 'w')
fileunit.puts 'OLR' + ': ' + gphys2.to_s
fileunit.puts 'OSR' + ': ' + gphys5.to_s
fileunit.close