# -*- coding: euc-jp -*-
require "numru/ggraph"
include NumRu

casenum = 27

gplw = GPhys::IO.open("DTempDtRadL.nc","DTempDtRadL")
gpsw = GPhys::IO.open("DTempDtRadS.nc","DTempDtRadS")

gplw = gplw.cut('lon'=>0).cut('lat'=>0).cut('time'=>0).copy
gplw = gplw * 60.0 * 60.0 * 24.0
gplw.long_name = 'Radiative heating rate'
gplw.units = 'K (day)-1'

gpsw = gpsw.cut('lon'=>0).cut('lat'=>0).cut('time'=>0).copy
gpsw = gpsw * 60.0 * 60.0 * 24.0
gpsw.long_name = 'Radiative heating rate'
gpsw.units = 'K (day)-1'


ps = 1013.00
p 'Surface pressure is assumed to be ' + ps.to_s + ' hPa'

plev = gplw.coord('sig').val
plev = plev * ps
varrayplev = VArray.new( plev,
                         { "long_name"=>'Pressure',
                           "units"=>'hPa' }, "plev" )
gplw.axis('sig').set_pos(varrayplev)
gpsw.axis('sig').set_pos(varrayplev)

gplwlblrtm = GPhys::IO.open("lblrtm/Case"+casenum.to_s+"hr.nc","hr")
z = gplwlblrtm.axis('level').pos.convert_units( Units['hPa'] )
gplwlblrtm.axis('level').set_pos(z)
gplwlblrtm = gplwlblrtm * 86400.0

iws= 4
DCL.gropn(iws)                            # iws : DCL出力デバイス．1,4:画面, 2:PS
DCL.glpset('lmiss',true)                  # DCLの欠損値処理を on に
DCL.sgpset('lclip',true)
GGraph.set_fig "window"=>[-5,5,ps,0], "viewport"=>[0.15,0.85,0.15,0.85]

GGraph.line gplw, true, 'exchange'=>true, 
  "annot"=>false,  # 右脇のアノテーション非表示
  "title"=>"Radiative heating rate", 'type'=>1, 'index'=>20, 
  'legend_vx'=>0.2, "legend"=>'Longwave'
GGraph.line gpsw, false, 'exchange'=>true, 
  "annot"=>false,  # 右脇のアノテーション非表示
  "title"=>"Radiative heating rate", 'type'=>1, 'index'=>40, 
  'legend_vx'=>0.2, "legend"=>'Shortwave'

GGraph.line gplwlblrtm, false, 'exchange'=>true, 
  "annot"=>false,  # 右脇のアノテーション非表示
  "title"=>"Radiative heating rate", 'type'=>2, 'index'=>20
#, 
#  'legend_vx'=>0.2, "legend"=>'Longwave (LBLRTM)'

DCL.grcls