# 月平均した降水量の緯度経度分布を描画するスクリプト
# dcpam と NCEP を 2 枚並べて描画する用
#
# * 2011/11/29 viewport などの細かい設定を修正
# * 2011/11/28 drawxy_clim_mon.rb をもとに作成

require "numru/ggraph"
include NumRu

######################################

def clim_mean( iyrs, iyre, imon, gphys, data )

  if data == "dcpam" then

     dataNumOfDay = 4
     daysOfMonth  = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]

     daysofyear = 0
     daysOfMonth.each do |i|
       daysofyear += i
     end

     time = gphys.coord('time').val

     for iyr in iyrs..iyre
       t = daysofyear*(iyr-1)
       for i in 1..imon-1
         t = t + daysOfMonth[i-1]
       end

       ts = t
       te = t + daysOfMonth[imon-1]

       ts = ts * dataNumOfDay
       te = te * dataNumOfDay - 1

       if iyr == iyrs
         gphystmpD = gphys.cut('time'=>time[ts]..time[te]).mean('time')
       else
         gphystmpD = gphystmpD + gphys.cut('time'=>time[ts]..time[te]).mean('time')
       end
         gphystmpD = gphystmpD / ( iyre - iyrs + 1 )
     end
 
     return gphystmpD

  elsif data == "NCEP" then

    hoursOfDay = 24
    daysofyear = 0

    time = gphys.coord('time').val

    for iyr in iyrs..iyre
 
     t  = (iyr - 1948) * 12
     t  = t + (imon - 1)
     ts = t
 
#      p iyrs, imon, ts

      if iyr == iyrs
        gphystmpN = gphys.cut('time'=>time[ts])
      else
        gphystmpN = gphystmpN + gphys.cut('time'=>time[ts])
      end
    end

     gphystmpN = gphystmpN / ( iyre - iyrs + 1 )

     return gphystmpN

  end

end

######################################

iyrsD = ARGV[4].to_i
iyreD = ARGV[5].to_i
iyrsN = ARGV[6].to_i
iyreN = ARGV[7].to_i

imon = ARGV[8].to_i

title = ARGV[9]

dirD    = ARGV[0]
dirN  = ARGV[1] 

#vname1 = 'SurfTemp'
#vname2 = 'Rain'
#vname3 = 'SoilMoist'
#vname4 = 'SurfSnow'

# 変数の設定

vnameD = ARGV[2]
vnameN = ARGV[3]

# gphys の準備

gphysD = GPhys::NetCDF_IO.open(dirD+'/'+vnameD+".nc", vnameD) 

tmD = gphysD.shape[2] # number of elements for time      dimension
timeD = gphysD.coord('time').val

gphysN = GPhys::NetCDF_IO.open(dirN+'/'+vnameN+ ".mon.mean"+".nc", vnameN) 
tmN = gphysN.shape[0] # number of elements for time      dimension
timeN = gphysN.coord('time').val



#< DCLのオープンと設定 >
DCL.gropn(2)
DCL.sgpset('lcntl', false)   # 制御文字を解釈しない
DCL.sgpset('lfull',true)     # 全画面表示
DCL.sgpset('lcorner', false)   # コーナーマークを書かない
DCL.uzfact(0.5)             # 座標軸の文字列サイズを 0.5 倍
DCL.sgpset('lfprop',true)    # プロポーショナルフォントを使う
DCL.udpset('lmsg',false)     # コンター間隔非表示

#< GGraph による 描画 >
vpt = NArray[0.2, 0.7, 0.01, 0.21]
vpt1 = (vpt + ([0.050]*2 + [0.32]*2) ).to_a
vpt2 = (vpt + ([0.050]*2 + [0.10]*2) ).to_a
GGraph.set_axes('ylabelint'=>30) # x 軸に値を打つ間隔

GGraph.set_fig 'itr'=>10 
GGraph.set_map 'coast_world'=>true, 'grid'=>false

data = "dcpam"
gphystmpD = clim_mean( iyrsD, iyreD, imon, gphysD, data )

data = "NCEP"
gphystmpN = clim_mean( iyrsN, iyreN, imon, gphysN, data )


if vnameD == 'SurfTemp'

  GGraph.tone( gphystmp, true, 
             'lev'=>[210,220,230,240,250,260,270,280,290,300],
             # レベル＆パターンを陽に指定
             'pat'=>[20999,30999,40999,45999,50999,55999,60999,70999,80999,85999,90999], 
             # パタンの方が1つ多→±∞まで
             'map_axes'=>true )
  #GGraph.contour( gphystmp, false, 'lev'=>[200,210,220,230,240,250,260,270,280,290,300], 'index'=>3 )      # 参考まで

elsif vnameD == 'Rain'

  # 1 枚目の描画の設定
  GGraph.set_fig('viewport'=>vpt1)
  GGraph.set_axes('xunits'=>'', 'yunits'=>'', 'xtitle'=>'')
  DCL.uzpset('labelxb',false)

  # 描画
  gphystmpD = gphystmpD / 2.5e6
  GGraph.tone( gphystmpD, true, 
             'lev'=>[0,1e-5,2e-5,4e-5,8e-5,12e-5,16e-5,20e-5,24e-5,28e-5,32e-5],
             # レベル＆パターンを陽に指定
             'pat'=>[1,20999,30999,35999,40999,50999,60999,65999,70999,75999,80999,85999], 
             # パタンの方が1つ多→±∞まで
             'map_axes'=>true,
             'annot'=>false, 'titl'=>'' )

  # 図の左上のタイトル
  DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','dcpam5',-1) ; DCL.uzpset('pad1',0.7)


  # 2 枚目の描画の設定
  GGraph.set_fig('viewport'=>vpt2,'new_frame'=>false)
  GGraph.set_axes('xunits'=>nil, 'yunits'=>'', 'xtitle'=>nil)
  DCL.uzpset('labelxb',true)

  # 描画
  GGraph.tone( gphystmpN, true, 
             'lev'=>[0,1e-5,2e-5,4e-5,8e-5,12e-5,16e-5,20e-5,24e-5,28e-5,32e-5],
             # レベル＆パターンを陽に指定
             'pat'=>[1,20999,30999,35999,40999,50999,60999,65999,70999,75999,80999,85999], 
             # パタンの方が1つ多→±∞まで
             'map_axes'=>true,
             'annot'=>false, 'titl'=>'' )

  # 図の左上のタイトル
  DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','NCEP',-1) ; DCL.uzpset('pad1',0.7)

  #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 )      # 参考まで

elsif vnameD == 'SoilMoist'

  GGraph.tone( gphystmp, true, 
             'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150],
             # レベル＆パターンを陽に指定
             'pat'=>[15999,20999,25999,30999,35999,40999,45999,50999,55999,60999,65999,70999,75999,80999,85999,90999,95999],
             # パタンの方が1つ多→±∞まで
             'map_axes'=>true )
  #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 )      # 参考まで

elsif vnameD == 'SurfSnow'

  GGraph.tone( gphystmp, true, 
             'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150],
             # レベル＆パターンを陽に指定
             'pat'=>[15999,20999,25999,30999,35999,40999,45999,50999,55999,60999,65999,70999,75999,80999,85999,90999,95999], 
             # パタンの方が1つ多→±∞まで
             'map_axes'=>true )
  #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 )      # 参考まで

elsif vnameD == 'OLRA'

  GGraph.tone( gphystmp, true, 
             'lev'=>[100,120,140,160,180,200,220,240,260,280,300],
             # レベル＆パターンを陽に指定
             'pat'=>[10999,15999,20999,25999,30999,35999,40999,50999,60999,70999,80999,90999],
             # パタンの方が1つ多→±∞まで
             'map_axes'=>true )
  #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 )      # 参考まで

end


# カラーバーの設定
GGraph.color_bar( 'voff'=>0.05, 'vcent'=>0.3,'units'=>'(kg/m/s)' )  

# 図のタイトルの設定
DCL::sgtxzv(0.5,vpt1[3]+0.06, 'Monthly mean precipitation '+'('+title+')',
            1.15*DCL.uzpget('rsizec2'),0,0,3)


DCL.grcls

