#!/usr/bin/env ruby
=begin
= 放射領域平均値時系列作図スクリプト(赤道, 中緯度)

   (e)放射収支(地表面正味長波フラックス)

=end

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

## 描画メソッドメイン

def plot_line_main(gphys_array)
  
  line_index_ary = Array.new; name_ary = Array.new

  default_index = 12
  
  line_hash = { "index"=> default_index }
  line_index_ary.push(default_index)
  name_ary.push(gphys_array[0].data.get_att("long_name"))
  
  x = gphys_array[0].coord(0).val
  y = gphys_array[0].data.val
  DCL.uulinz(x, y, 1, default_index)
  gphys_array.size.times{ |num|
      unless num == 0
	line_index = (num*10 + default_index)
	line_index_ary.push(num*10 +default_index)
	name_ary.push(gphys_array[num].data.get_att("long_name"))
	x = gphys_array[num].coord(0).val
	y = gphys_array[num].data.val
	DCL.uulinz(x, y, 1, line_index)
      end
  }
  
  return name_ary, line_index_ary
  
end



##-----------------------
#  ラインインデックスの種類を表示

def __show_line_index(str_ary,line_index_ary, x=0.15, y=0.12) 
  charsize = 1.0 * DCL.uzpget('rsizec1')
  dvx = 0.01
  dvy = charsize*1.5
  raise TypeError,"Array expected" if ! str_ary.is_a?(Array)
  vxmin,vxmax,vymin,vymax = DCL.sgqvpt
  vx = x
  vy = y + charsize/2
  str_ary.size.times{|num|
    DCL::sgtxzv(vx,vy,"--- #{str_ary[num]}", 
		  charsize, 0, -1, line_index_ary[num])
      vy -= dvy
      if num == 2
	vx = 0.30
	vy = 0.12 - charsize/2
      end
  }
    nil
  end



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

# 大気上端上向き長波放射
nlwrs_long_name = "NLR Surface"
nlwrs_nh_law_00_10_long_name = nlwrs_long_name + " ( EQ - N10)"
nlwrs_nh_law_10_20_long_name = nlwrs_long_name + " (N10 - N20)"
nlwrs_nh_law_20_30_long_name = nlwrs_long_name + " (N20 - N30)"
nlwrs_sh_law_00_10_long_name = nlwrs_long_name + " ( EQ - S10)"
nlwrs_sh_law_10_20_long_name = nlwrs_long_name + " (S10 - S20)"
nlwrs_sh_law_20_30_long_name = nlwrs_long_name + " (S20 - S30)"
nlwrs_nh_midlat_long_name = nlwrs_long_name + " (N30 - N60)"
nlwrs_sh_midlat_long_name = nlwrs_long_name + " (S30 - S60)"


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

##  ガウス重み
path_gw = "../../../GAUSSIAN_WEIGHT/gaussian_weight.nc"
gp_gw   = GPhys::NetCDF_IO.open(path_gw, 'gw')

axmonth = []
i = 0

nlwrs_name = ""; nlwrs_unit = ""

nlwrs_nh_law_00_10 = []
nlwrs_nh_law_10_20 = []
nlwrs_nh_law_20_30 = []
nlwrs_sh_law_00_10 = []
nlwrs_sh_law_10_20 = []
nlwrs_sh_law_20_30 = []
nlwrs_nh_mid = []
nlwrs_sh_mid = []


year[0].upto(year[1]) do |y|
  month.each do |m|

    # NetCDFVar => GPhys 
    path_nlwrs = "../../../NLWRS.NCEP/NLWRS.#{y}.NCEP/NLWRS_#{y}-#{m}_NCEP.nc" # 地表面正味短波放射

    gp_nlwrs = GPhys::NetCDF_IO.open(path_nlwrs, 'nlwrs').mean('lon')

    # 時系列 VArray 作り
    ## nlwrs
    nlwrs_nh_law_00_10 << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=>  0..10  ) ).sum.val )
    nlwrs_nh_law_10_20 << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=> 10..20  ) ).sum.val )
    nlwrs_nh_law_20_30 << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=> 20..30  ) ).sum.val )
    nlwrs_sh_law_00_10 << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=>  0..-10 ) ).sum.val )
    nlwrs_sh_law_10_20 << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=>-10..-20 ) ).sum.val )
    nlwrs_sh_law_20_30 << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=>-20..-30 ) ).sum.val )

    nlwrs_nh_mid << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val )
    nlwrs_sh_mid << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val )

    nlwrs_name = gp_nlwrs.data.name
    nlwrs_unit = gp_nlwrs.data.units.to_s

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


### for intensity

  # 日付軸作成
  namonth = NArray.to_na(axmonth)                 # 日付軸の NArray
  vamonth = VArray.new(namonth).rename!("month")  # VArray 化. 名前は "month"
  vamonth.set_att('long_name','month')            # 属性設定. long_name と
  vamonth.set_att('units','month since 1979 JAN') #           units     を設定.
  grid = Grid.new( Axis.new(true).set_cell_guess_bounds(vamonth).set_pos_to_center )

  # nlwrs を GPhys 化
  ## VArray -> GPhys
  gp_nlwrs_nh_law_00_10 = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_nh_law_00_10)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=>  0..10  ).sum.val 
  gp_nlwrs_nh_law_10_20 = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_nh_law_10_20)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=> 10..20  ).sum.val 
  gp_nlwrs_nh_law_20_30 = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_nh_law_20_30)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=> 20..30  ).sum.val 
  gp_nlwrs_sh_law_00_10 = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_sh_law_00_10)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=>  0..-10 ).sum.val 
  gp_nlwrs_sh_law_10_20 = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_sh_law_10_20)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=>-10..-20 ).sum.val 
  gp_nlwrs_sh_law_20_30 = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_sh_law_20_30)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=>-20..-30 ).sum.val 

  gp_nlwrs_nh_mid       = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_nh_mid)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=> 60..30  ).sum.val 
  gp_nlwrs_sh_mid       = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_sh_mid)).rename!(nlwrs_name)) / gp_gw.cut( 'lat'=>-30..-60 ).sum.val 
  ## データの属性設定
  gp_nlwrs_nh_law_00_10.data.set_att('long_name', nlwrs_nh_law_00_10_long_name)
  gp_nlwrs_nh_law_10_20.data.set_att('long_name', nlwrs_nh_law_10_20_long_name)
  gp_nlwrs_nh_law_20_30.data.set_att('long_name', nlwrs_nh_law_20_30_long_name)
  gp_nlwrs_sh_law_00_10.data.set_att('long_name', nlwrs_sh_law_00_10_long_name)
  gp_nlwrs_sh_law_10_20.data.set_att('long_name', nlwrs_sh_law_10_20_long_name)
  gp_nlwrs_sh_law_20_30.data.set_att('long_name', nlwrs_sh_law_20_30_long_name)
  gp_nlwrs_nh_mid.data.set_att('long_name', nlwrs_nh_midlat_long_name)
  gp_nlwrs_sh_mid.data.set_att('long_name', nlwrs_sh_midlat_long_name)
  gp_nlwrs_nh_law_00_10.data.set_att('units', nlwrs_unit)
  gp_nlwrs_nh_law_10_20.data.set_att('units', nlwrs_unit)
  gp_nlwrs_nh_law_20_30.data.set_att('units', nlwrs_unit)
  gp_nlwrs_sh_law_00_10.data.set_att('units', nlwrs_unit)
  gp_nlwrs_sh_law_10_20.data.set_att('units', nlwrs_unit)
  gp_nlwrs_sh_law_20_30.data.set_att('units', nlwrs_unit)
  gp_nlwrs_nh_mid.data.set_att('units', nlwrs_unit)
  gp_nlwrs_sh_mid.data.set_att('units', nlwrs_unit)

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

##
# 事前準備

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

##
# お絵描きメイン

# 縦軸の最大描画範囲(W/m^2)
max = +100
min = 0
#max = 135
#min = 125

# ビューポート設定
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)            # 日日数
#GGraph.set_axes("date?"=>["x", idate, nd])        
                                          # 日付軸を有効にする


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 軸の単位の位置を下方へ 


### 1 枚目 : 北半球低緯度時系列

DCL::grfrm                   # ページ確定 

DCL::grswnd(0.0, nd, min, max)
DCL::grsvpt(*vpt)
DCL::grstrn(1)
DCL::grstrf

DCL::ucxacl('B', idate, nd)
DCL::ucxacl('T', idate, nd)

DCL::uyaxdv('L', max/10, max/5)
DCL::uyaxdv('R', max/10, max/5)

DCL::uxmttl('T', 'Mean net long wave radiation flux at NH lawlat.(1979-2003)', 0.0)
DCL::uysttl('L', "radiation(#{nlwrs_unit.to_s})", 0.0)
DCL::uxsttl('B', 'YEAR', 0.0)
# DCL::uscset('cyunit', nlwrs_unit.to_s ) # y 軸の単位

names, idxs = plot_line_main([ 
			       gp_nlwrs_nh_law_00_10,
			       gp_nlwrs_nh_law_10_20,
			       gp_nlwrs_nh_law_20_30
                              ])

__show_line_index(names,idxs)

### 2 枚目 : 南半球低緯度時系列

DCL::grfrm                   # ページ確定 

DCL::grswnd(0.0, nd, min, max)
DCL::grsvpt(*vpt)
DCL::grstrn(1)
DCL::grstrf

DCL::ucxacl('B', idate, nd)
DCL::ucxacl('T', idate, nd)

DCL::uyaxdv('L', max/10, max/5)
DCL::uyaxdv('R', max/10, max/5)

DCL::uxmttl('T', 'Mean net long wave radiation flux at SH lawlat.(1979-2003)', 0.0)
DCL::uysttl('L', "radiation(#{nlwrs_unit.to_s})", 0.0)
DCL::uxsttl('B', 'YEAR', 0.0)
# DCL::uscset('cyunit', nlwrs_unit.to_s ) # y 軸の単位

names, idxs = plot_line_main([ 
			       gp_nlwrs_sh_law_00_10,
			       gp_nlwrs_sh_law_10_20,
			       gp_nlwrs_sh_law_20_30
                              ])

__show_line_index(names,idxs)

### 3 枚目 : 中緯度時系列

DCL::grfrm                   # ページ確定 

DCL::grswnd(0.0, nd, min, max)
DCL::grsvpt(*vpt)
DCL::grstrn(1)
DCL::grstrf

DCL::ucxacl('B', idate, nd)
DCL::ucxacl('T', idate, nd)

DCL::uyaxdv('L', max/10, max/5)
DCL::uyaxdv('R', max/10, max/5)

DCL::uxmttl('T', 'Mean net long wave radiation flux at midlat.(1979-2003)', 0.0)
DCL::uysttl('L', "radiation(#{nlwrs_unit.to_s})", 0.0)
DCL::uxsttl('B', 'YEAR', 0.0)
# DCL::uscset('cyunit', nlwrs_unit.to_s ) # y 軸の単位

names, idxs = plot_line_main([ 
			       gp_nlwrs_nh_mid, gp_nlwrs_sh_mid
                              ])

__show_line_index(names,idxs)


DCL.grcls


##################################################
# 画像ファイル名を変更
File.rename('dcl_001.ps', 'NLWRS_NH_LAW_1979-2003_TIME-SERIES_NCEP.ps')
File.rename('dcl_002.ps', 'NLWRS_SH_LAW_1979-2003_TIME-SERIES_NCEP.ps')
File.rename('dcl_003.ps', 'NLWRS_MID_1979-2003_TIME-SERIES_NCEP.ps')