#!/usr/bin/env ruby
=begin

=draw_energy.rb

  * エネルギーの時間変化の折れ線図を描く
  * 初期のエネルギーに対する比を描いている

=end

def energyaccml(gp_tend)

timestep = 600.0
time = gp_tend.axis(0)
len = gp_tend.data.length
ary  = NArray.sfloat(len).fill(1.0)

i = 0
for tend in gp_tend.data.val
  unless i==0
    ary[i] = ary[i-1] + tend * (time.pos.val[i] - time.pos.val[i-1])
  else
    ary[i] = tend * time.pos.val[i]
  end 
  i=i+1
end
data = VArray.new( ary,
                   {"long_name"=>"energy change by each components", "units"=>"J"},
                   "Energy" )
gp_accml = GPhys.new(Grid.new(time), data) 
return gp_accml

end

require "numru/ggraph"
require "numru/dcl"

include NumRu

Timestep = 600.0 # 出力時間間隔 [sec]

#file = ARGV[0]
file = ["nc/arare_diag-odaka1998-dx100-20051130.nc",
        "nc/arare_diag-odaka1998-dx100-20051130-2.nc"]

# open file and initialize GPhys object

gp_rad = GPhys::IO.open(file, "EnergyTendRadTotal")
gp_dis = GPhys::IO.open(file, "EnergyTendDisTotal")
gp_sfh = GPhys::IO.open(file, "EnergyTendSfcHeatFluxTotal")
gp_sfm = GPhys::IO.open(file, "EnergyTendSfcMomentumFluxTotal")
gp_ncs = GPhys::IO.open(file, "EnergyTendNonconsvTotal")
gp_tbh = GPhys::IO.open(file, "EnergyTendPotTempTurbTotal")
gp_tbm = GPhys::IO.open(file, "EnergyTendVelTurbTotal")
gp_bsf = GPhys::IO.open(file, "EnergyTendPotTempBasicTotal")

# create new GPhys objects

gp_rad2 = energyaccml(gp_rad)
gp_dis2 = energyaccml(gp_dis)
gp_sfh2 = energyaccml(gp_sfh)
gp_sfm2 = energyaccml(gp_sfm)
gp_ncs2 = energyaccml(gp_ncs)
gp_tbh2 = energyaccml(gp_tbh)
gp_tbm2 = energyaccml(gp_tbm)
gp_bsf2 = energyaccml(gp_bsf)

gp_accml_total = gp_rad2 + gp_dis2 + gp_sfh2 + gp_sfm2  \
                 + gp_ncs2 + gp_tbh2 + gp_tbm2 + gp_bsf2

# draw figure

DCL.gropn(4)
DCL.sgpset("lcntl", false)
DCL.uzfact(0.7)
GGraph.line(gp_accml_total, true, "index"=>13, "max"=>6.0e+9, "min"=>-1.0e+9)
GGraph.line(gp_bsf2, false, "index"=>23 )
GGraph.line(gp_rad2, false, "index"=>33)
GGraph.line(gp_dis2, false, "index"=>63 )
GGraph.line(gp_sfh2, false, "index"=>43 )
GGraph.line(gp_sfm2, false, "index"=>43, "type"=>2 )
GGraph.line(gp_tbh2, false, "index"=>73 )
GGraph.line(gp_tbm2, false, "index"=>73, "type"=>2 )
GGraph.line(gp_ncs2, false, "index"=>83 )
DCL.grcls