#!/usr/bin/env ruby
=begin
= mkfig-CF.rb -- 地球大気における放射強制力の緯度分布
=end

##########  設定部分 ここから #################
year = [1980, 2000]
range = [nil, nil]

month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]   # 初期値 
lat   = [-30, 30]

var = "prate"
var_dir = "PRATE"

# ps に落す 
$OPT_f = true

# 絵 タイトル
fig_title      = "Precipitation at Equatorial Surface"

# line タイトル
line_title_CMAP = "CMAP Precipitation at Equatorial Surface"
line_title_NCEP = "NCEP Precipitation at Equatorial Surface"
line_title_CMAP_MODEL = "CMAP-MODEL Precipitation at Equatorial Surface"

# データソース
src_CMAP = "CMAP"
src_NCEP = "NCEP"

# y 軸タイトル
y_axis_title = "Precipitation rate"

vp = [0.15, 0.85, 0.2, 0.55]


if ARGV == 4
  year = ARGV[0..1]
  range = ARGV[2..3]
end

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

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

module NumRu
  class NetCDFVar
    alias get get_with_miss_and_scaling
    def self.put(*args); put_with_miss_and_scaling(*args) end
  end
end


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




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

def make_gphys(var_dir, var, lat, year, month, figtitle, yaxtitle, line_title, src)

  axmonth = []
  i = 0
  gprate = []
  va_name = ""
  va_unit = ""
  
  year[0].upto(year[1]) do |y|
    month.each do |m|
      if src == "CMAP"
	path = "../../#{src}/#{var_dir}.#{src}/"
      else
	path = "../../#{src}/#{var_dir}.#{src}/#{var_dir}.#{y}.#{src}/"
      end
      name = "#{var_dir}_#{y}-#{m}_#{src}.nc"
      
      gphys = GPhys::NetCDF_IO.open(path+name, var)                        
      
      gphys = gphys.cut("lat"=>lat[0]..lat[1])
      
      gprate << gphys.mean(0).mean
      
      va_name = gphys.data.get_att("long_name")
      va_unit = gphys.data.get_att("units")
      
      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 = GPhys.new(grid, VArray.new(NArray.to_na(gprate)).rename!(yaxtitle).set_att("line_name", line_title) )
  
  ## データの属性設定
  gphys.data.set_att('units',va_unit)
  
  # 軸の属性設定
  gphys.coord(0).set_att('long_name','month')
  gphys.coord(0).set_att('units','month since 1980 JAN')

  return gphys
  
end

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


gphys_cmap = make_gphys(var_dir, var, lat, year, month, fig_title, y_axis_title, line_title_CMAP, src_CMAP)
gphys_ncep = make_gphys(var_dir, var, lat, year, month, fig_title, y_axis_title, line_title_NCEP, src_NCEP)
gphys_cmap_model = make_gphys(var_dir, var, lat, year, month, fig_title, y_axis_title, line_title_CMAP_MODEL, src_CMAP)


##
# 事前準備

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_cmap, gphys_ncep, gphys_cmap_model], "min"=>range[0], "max"=>range[1])

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

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