#!/usr/bin/env ruby
=begin
= mkfig-PR_EQ.rb -- 赤道域降雨量の時系列図
=end

##########  設定部分 ここから #################
# 対象とする年範囲
year = [1980, 2001]
# 対象とする月
month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]   # 初期値 
# netCDF 変数名
var = "prate"
# ディレクトリ名
var_dir = "PRATE"
# 対象緯度範囲
lat_range_north = [60,30]
lat_range_south = [-30,-60]
# 表示する単位
va_unit = "kg/s"

# ps に落す 
$OPT_f = true

# タイトル
fig_title = "Precipitation at Midlatitude Surface"
fig_title_north = "Precipitation at Midlatitoude Surface on northsphere"
fig_title_south = "Precipitation at Midlatitude Surface on southsphere"
# y 軸タイトル
y_axis_title = "Precipitation"

vp = [0.15, 0.85, 0.2, 0.55]

# 地球半径(m)
R = 6378e3

##########  設定部分 ここまで #############

require "getopts"        # for option_parse
require "numru/netcdf"
require "numru/ggraph"
include NumRu
require "libdraw-n"
include Draw
require "libgphys-n"


# set User Path for dcldatabase
DCL.glcset('DUPATH','/home/daktu32/.dcldir/')     


########## サブルーチン ここから ##############

##
# 各格子点に対応する面積要素の面積を求めるメソッド
# m^-2 で規格化されている降雨量を修正する.

def cal_areas(r, lat_deg, lon_d_deg, lat_d_deg)
 a = lat_deg/360.0*2*Math::PI
 b = lon_d_deg/360.0*2*Math::PI
 c = lat_d_deg/360.0*2*Math::PI
 area = r*r*b*(Math::sin(a+c/2)-Math::sin(a-c/2))
 return area
end


########## サブルーチン ここまで ############## 


################################################################
#                        make gphys 
################################################################


## make compisit-GPhys object 
# make array which set the gphys objects.
axmonth = []
i = 0
gprate_north = []
gprate_south = []

year[0].upto(year[1]) do |y|
  month.each do |m|
    path = "../../#{var_dir}.NCEP/#{var_dir}.#{y}.NCEP/"
    name = "#{var_dir}_#{y}-#{m}_NCEP.nc"

    gphys = GPhys::NetCDF_IO.open(path+name, var)                        
    lon_d_deg = (gphys.coord(0).val[0]-gphys.coord(0).val[1]).abs
    lat_d_deg = (gphys.coord(1).val[0]-gphys.coord(1).val[1]).abs

    int_north = []
    int_south = []

    gphys_north = gphys.cut("lat"=>lat_range_north[0]..lat_range_north[-1])
    gphys_south = gphys.cut("lat"=>lat_range_south[0]..lat_range_south[-1])

    gphys_north.coord(1).val.each {|lat|
    darea = cal_areas(R, lat, lon_d_deg, lat_d_deg)
    int_north << gphys_north.cut("lat"=>lat).data.val.sum*darea  # 緯度円に沿って積分し, グリッドに対応する面積を掛ける
    }

    gphys_south.coord(1).val.each {|lat|
    darea = cal_areas(R, lat, lon_d_deg, lat_d_deg)
    int_south << gphys_south.cut("lat"=>lat).data.val.sum*darea  # 緯度円に沿って積分し, グリッドに対応する面積を掛ける
    }
    
    total_int_north = NArray.to_na(int_north).sum
    total_int_south = NArray.to_na(int_south).sum

    gprate_north << total_int_north
    gprate_south << total_int_south

    axmonth << DCL::dateg3(year[0],1,1,y.to_i,m.to_i,1)
    i += 1
  end
end

axmonth = NArray.to_na(axmonth)
axmonth = Axis.new(true).set_cell_guess_bounds(VArray.new(axmonth).rename!("month")).set_pos_to_center

grid = Grid.new(axmonth)

gphys_north = GPhys.new(grid, VArray.new(NArray.to_na(gprate_north)).rename!(y_axis_title).set_att("line_name", fig_title_north) )
gphys_south = GPhys.new(grid, VArray.new(NArray.to_na(gprate_south)).rename!(y_axis_title).set_att("line_name", fig_title_south) )

## データの属性設定
gphys_north.data.set_att('units',va_unit)
gphys_south.data.set_att('units',va_unit)

# 軸の属性設定
gphys_north.coord(0).set_att('long_name','month')
gphys_south.coord(0).set_att('units','month since 1980 JAN')


################################################################
#                        描画ルーチン
################################################################

##
# 事前準備

DCL.uscset('cyspos', 'B' )              # y 軸の単位の位置を下方へ 
rsizel2 = DCL.uzrget('rsizel2')         # 現在のラベルサイズを取得
DCL.uzrset('rsizel2', rsizel2*0.42 )    # ラベルサイズをデフォルトの 0.5 倍に

##
# お絵描きメイン

Draw::mkwin(vp)                         # ウィンドウオープン

idate = (year[0].to_s+"0101").to_i      # 日付 0 日. 
ndate = (year[-1].to_s+"1201").to_i     # 日付 最後の日. 
nd = DCL.dateg1(idate,ndate)            # 日日数
GGraph.set_axes("date?"=>["x", idate, nd])        
                                        # 日付軸を有効にする
Draw::plot_line_main([gphys_south, gphys_north])

DCL::uxsttl("b", "YEAR", 0)             # x 軸タイトルを "YEAR" に 
DCL::uxmttl("t", fig_title, -1)         # タイトル

Draw::clwin                             # ウィンドウクローズ