#!/usr/bin/env ruby
=begin
= mkfig_25time_series.rb -- make time-series figure of zonal wind  at 200 hPa

=end

require "getopts"        # for option_parse
require "numru/ggraph"
include NumRu
include Misc::EMath


################################################################
#                        make gphys 
################################################################
year = [1979, 2003]
month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]   # 初期値 

lat_range = 50..-50 # 緯度範囲

## make compisit-GPhys object 
# make array which set the gphys objects.

axmonth = []
i = 0

uwnd_name = ""  # 下記のイテレータ中で変数名を取得する. スコープの関係上ここで宣言.
uwnd_unit = ""  # 

uwnd_at_200hPa = []

year[0].upto(year[1]) do |y|
  month.each do |m|
    path = "../../../UWND.NCEP/UWND.#{y}.NCEP/UWND_#{y}-#{m}_NCEP.nc" # zonal wind

    gp_uwnd = GPhys::NetCDF_IO.open(path, 'uwnd').cut('lat'=>lat_range)

    va_uwnd_at_200hPa = gp_uwnd.cut( 'level'=>200  ).mean('lon').data.copy
                     
    uwnd_at_200hPa << va_uwnd_at_200hPa.reshape!(1, *va_uwnd_at_200hPa.shape)
		                           
    uwnd_name =      gp_uwnd.data.name
    uwnd_unit =      gp_uwnd.data.get_att("units")

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


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

    
### for width

  va_uwnd_at_200hPa = VArrayComposite.new( NArray.to_na( uwnd_at_200hPa ) )
  va_uwnd_at_200hPa.set_att('long_name', 'Deviation of zonal wind speed at 200 hPa (1979-2003)')

  uwnd_path = "../../25MEANS.NCEP/25MEANS.1979-2003.NCEP/100hPa_PHYS_1979-2003_ANNUAL_NCEP.nc"
  gp_uwnd_at_200hPa_climat = GPhys::IO.open(uwnd_path, "uwnd").cut('level'=>200, 'lat'=>lat_range)
  axlat = gp_uwnd_at_200hPa_climat.axis('lat')
  gp_uwnd_at_200hPa = GPhys.new( Grid.new(axmonth, axlat), va_uwnd_at_200hPa )
  gp_uwnd_deviation = gp_uwnd_at_200hPa - gp_uwnd_at_200hPa_climat

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

##
# 事前準備

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

##
# お絵描きメイン

# ビューポート設定
vpt = NArray[0.15,0.85,0.23,0.58]

p idate = (year[0].to_s+"0101").to_i      # 日付 0 日. 19790101
ndate = (year[-1].to_s+"1201").to_i       # 日付 最後の日. 20031231
p nd = DCL.dateg1(idate,ndate)            # 日日数


DCL.swpset('lsep',  true)    # ページ毎別々のファイルに落す
DCL.gropn(2)
DCL.sgpset('lcntl', false)   # 
DCL.sgpset('lfull',true)     # フルスクリーン
DCL.uzfact(0.55)             # 
DCL.sgpset('lcorner',false)  # コーナーを取っ払う 
DCL.sgpset('lfprop',true)    # 
DCL.udpset('lmsg',false)     # 
DCL.uscset('cyspos', 'B' )   # y 軸の単位の位置を下方へ 



## set variable param 

level   = NArray[-10, -5, 0, 5, 10]
pattern = NArray[30999, 35999, 40999, 70999, 75999, 85999]
cont_min = -10
cont_max = +10
cont_int = 0.5 

# get coord info
missval = gp_uwnd_at_200hPa_climat.data.get_att('missing_value')[0]
before = DCLExt.gl_set_params('lmiss'=>true,'rmiss'=>missval)

ymax = gp_uwnd_deviation.coord('lat').max
ymin = gp_uwnd_deviation.coord('lat').min

###  200 hPa 面時系列
  DCL::grfrm                   # ページ確定 

  DCL::grswnd(0.0, nd, ymin, ymax)
  DCL::grsvpt(*vpt)
  DCL::grstrn(1)
  DCL::grstrf

  ## トーン, コンター描画
  DCLExt::ue_set_tone(level, pattern)
  DCLExt::ud_set_linear_levs(gp_uwnd_deviation.data.val, 'int'=>cont_int, 'min'=>cont_min, 'cont_max'=>cont_max)
  DCL::uetone(gp_uwnd_deviation.data.val[true, -1..0])
  DCL::udcntz(gp_uwnd_deviation.data.val[true, -1..0])

  ## 軸描画

  ### 日付軸
  DCL::ucxacl('B', idate, nd)
  DCL::ucxacl('T', idate, nd)
  ### y 軸ラベル設定
  DCL::uyaxdv('L', 10, 20)
  DCL::uyaxdv('R', 10, 20)
  DCL::uxsttl('T', "(#{gp_uwnd_deviation.units.to_s})", 1)
  DCL::uxmttl('T',  gp_uwnd_deviation.data.get_att('long_name'), 0.0)
  DCL::uysttl('L', "Latitude", 0.0)
  DCL::uxsttl('B', 'YEAR', 0.0)

  ## color bar
  DCLExt::color_bar('constwidth'=>true)

  DCL::ueitlv

DCL.grcls

##################################################
# 画像ファイル名を変更
File.rename('dcl_001.ps', 'UWND_1973-2003_TIME-SERIES_200hPa_NCEP_DEVIATION.ps')