#!/usr/bin/env ruby
=begin

=draw_thetatend.rb

  * 診断量出力ファイルから温位の式の各項(以下参照)の寄与を描く
    * 移流項
    * 乱流拡散項(地表面フラックス含む, 数値拡散は含めない)
    * 放射冷却項
    * 散逸加熱項
    * 以上の合計

=end

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

include NumRu

Daytime = UNumeric.new(86400.0, "s/day")

file = "/home/kitamo/arare3/arare_diag-odaka1998-dx200-20050917.nc"  # NetCDF file
Starttime = 39600     
Endtime = 43200


# open file and initialize GPhys object

gp_adv = GPhys::IO.open(file, "PotTempTendAdvZ")
gp_src = GPhys::IO.open(file, "PotTempTendSrcZ")
gp_turb = GPhys::IO.open(file, "PotTempTendTurbZ")
gp_rad = GPhys::IO.open(file, "PotTempTendRadZ")
gp_dis = GPhys::IO.open(file, "PotTempTendDisZ")
gp_sfc = GPhys::IO.open(file, "PotTempTendSfcZ")

# unit convertion

gp_adv = gp_adv * Daytime
gp_src = gp_src * Daytime
gp_turb = gp_turb * Daytime
gp_rad = gp_rad * Daytime
gp_dis = gp_dis * Daytime
gp_sfc = gp_sfc * Daytime

# cut

gp_adv = gp_adv.cut('t'=>Starttime..Endtime)
gp_src = gp_src.cut('t'=>Starttime..Endtime)
gp_turb = gp_turb.cut('t'=>Starttime..Endtime)
gp_rad = gp_rad.cut('t'=>Starttime..Endtime)
gp_dis = gp_dis.cut('t'=>Starttime..Endtime)
gp_sfc = gp_sfc.cut('t'=>Starttime..Endtime)

# mean

gp_adv = gp_adv.mean('t')
gp_src = gp_src.mean('t')
gp_turb = gp_turb.mean('t')
gp_rad = gp_rad.mean('t')
gp_dis = gp_dis.mean('t')
gp_sfc = gp_sfc.mean('t')

# total tendency calculation

gp_diff = gp_turb + gp_sfc
#gp_adv =  gp_adv + gp_src
gp_total = gp_adv + gp_turb + gp_rad + gp_dis + gp_sfc + gp_src

# drawing

DCL.gropn(4)
DCL.sgpset('lcntl', false)
DCL.uzfact(0.7)
GGraph.line( gp_total, true, 'exchange'=>true, 'title'=>'POTENTIAL TEMPARETURE TENDENCY',
             'max'=>150, 'min'=>-150, 'index'=>12 )
GGraph.line( gp_adv, false, 'exchange'=>true, 'index'=>22 )
GGraph.line( gp_src, false, 'exchange'=>true, 'index'=>92 )
GGraph.line( gp_rad, false, 'exchange'=>true, 'index'=>32 )
#GGraph.line( gp_turb, false, 'exchange'=>true, 'index'=>42 )
GGraph.line( gp_dis, false, 'exchange'=>true, 'index'=>62 )
#GGraph.line( gp_sfc, false, 'exchange'=>true, 'index'=>72 )
GGraph.line( gp_diff, false, 'exchange'=>true, 'index'=>72 )
DCL.grcls