#!/usr/bin/env ruby

#  * Drawing tendency of potential temperature
#    * advection
#    * turbulent diffusion
#    * radiation
#    * dissipation
#    * total
#  * Usage
#    * ruby draw_thetatend.rb --t1 aaaa --t2 bbbb --range yyy:zzz
#    * aaaa is the time of beginning of averaging. 
#    * aaaa is the time of end of averaging. 
#  * This ruby script is partially based on the gpview script.
#


require "numru/ggraph" # ggraph library
require "numru/dcl"    # dcl library
require "getoptlong"   # A library for analysis of command line option
include NumRu


def __split_range(range)

  if /(.*):(.*)/ =~ range
    if $1 == ""
      min = nil
    else
      min = $1.to_f
    end
    if $2 == ""
      max = nil
    else
      max = $2.to_f
    end
  elsif range == nil
    min = max = nil
  else
    raise "invalid range: variable subset specification error. split
range with
':'\n\n"
  end

  return min, max
end

parser = GetoptLong.new
parser.set_options( ['--t1', GetoptLong::REQUIRED_ARGUMENT],
                    ['--t2', GetoptLong::REQUIRED_ARGUMENT],
                    ['--range', GetoptLong::REQUIRED_ARGUMENT]
                  )

begin
  parser.each_option do |name, arg|
    eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'"
  end
    rescue
  exit(1)
end

min_range,  max_range  = __split_range($OPT_range)


T1 = $OPT_t1.to_f # time
T2 = $OPT_t2.to_f # time
DelT = T2 - T1
Z1 = 27000
Z2 = 30000

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

file = "MarsCond_Zprof.nc"  # NetCDF file

# open file and initialize GPhys object

gp_adv = GPhys::IO.open(file, "DensCloudAdv")
gp_turb = GPhys::IO.open(file, "DensCloudTurb")
gp_cond = GPhys::IO.open(file, "DensCloudCond")
gp_diff = GPhys::IO.open(file, "DensCloudDiff")
gp_fill = GPhys::IO.open(file, "DensCloudFill")
gp_fill2 = GPhys::IO.open(file, "DensCloudFillZero")


# unit convertion

#gp_adv = gp_adv * Daytime
#gp_turb = gp_turb * Daytime
#gp_rad = gp_rad * Daytime
#gp_dis = gp_dis * Daytime
#gp_cond = gp_cond * Daytime

# cut

gp_adv = gp_adv.cut('t'=>T1..T2,'z'=>Z1..Z2)
gp_turb = gp_turb.cut('t'=>T1..T2,'z'=>Z1..Z2)
gp_cond = gp_cond.cut('t'=>T1..T2,'z'=>Z1..Z2)
gp_diff = gp_diff.cut('t'=>T1..T2,'z'=>Z1..Z2)
gp_fill = gp_fill.cut('t'=>T1..T2,'z'=>Z1..Z2)
gp_fill2 = gp_fill2.cut('t'=>T1..T2,'z'=>Z1..Z2)

#gp_adv = gp_adv.cut('t'=>T1..T2)
#gp_turb = gp_turb.cut('t'=>T1..T2)
#gp_cond = gp_cond.cut('t'=>T1..T2)
#gp_diff = gp_diff.cut('t'=>T1..T2)
#gp_fill = gp_fill.cut('t'=>T1..T2)
#gp_fill2 = gp_fill2.cut('t'=>T1..T2)

gp_adv = gp_adv.sum('t')
gp_turb = gp_turb.sum('t')
gp_cond = gp_cond.sum('t') 
gp_diff = gp_diff.sum('t') 
gp_fill = gp_fill.sum('t') 
gp_fill2 = gp_fill2.sum('t') 

# mean

#gp_adv = gp_adv.mean('t')
#gp_turb = gp_turb.mean('t')
#gp_cond = gp_cond.mean('t') 
#gp_diff = gp_diff.mean('t') 
#gp_fill = gp_fill.mean('t') 
#gp_fill2 = gp_fill2.mean('t') 

gp_adv = gp_adv.mean('x')
gp_turb = gp_turb.mean('x')
gp_cond = gp_cond.mean('x') 
gp_diff = gp_diff.mean('x') 
gp_fill = gp_fill.mean('x') 
gp_fill2 = gp_fill2.mean('x') 

# total tendency calculation

#gp_diff = gp_turb + gp_sfc
#gp_adv =  gp_adv + gp_src
gp_total = gp_adv + gp_turb + gp_cond + gp_diff + gp_fill + gp_fill2
gp_others = gp_turb + gp_fill + gp_fill2
#gp_others = gp_fill2
#gp_adv = gp_adv + gp_fill

#gp_adv = gp_adv / DelT
#gp_turb = gp_turb / DelT
#gp_cond = gp_cond / DelT
#gp_diff = gp_diff / DelT
#gp_fill = gp_fill / DelT
#gp_fill2 = gp_fill2 / DelT
#gp_total = gp_total / DelT
#gp_others = gp_others / DelT

gp_adv = gp_adv 
gp_turb = gp_turb 
gp_cond = gp_cond 
gp_diff = gp_diff 
gp_fill = gp_fill 
gp_fill2 = gp_fill2 
gp_total = gp_total 
gp_others = gp_others 

# drawing

DCL.gropn(4)
DCL.sgpset('lcntl', false)
DCL.uzfact(0.7)
DCL.sgpset('lfull', true)      # use full area in the window
DCL.sgpset('lfprop',true)      # use proportional font


GGraph::set_fig('view'=>[0.15,0.85,0.20,0.55])
GGraph.set_axes('xunits'=>'kg m^-3 s^-1 ','yunits'=>'m','xtitle'=>'','ytitle'=>'Z-coordinate')
#GGraph.line( gp_total, true, 'exchange'=>true, \
GGraph.line( gp_adv, true, 'exchange'=>true, \
'title'=>'CLOUD DENSITY TENDENCY', \
'min'=> __split_range($OPT_range)[0], \
'max'=> __split_range($OPT_range)[1], \
#'index'=>14 )  # black
'index'=>22 )  # red
#GGraph.set_axes('xunits'=>'kg m^-1 s^-1 ','yunits'=>'m','xtitle'=>'','ytitle'=>'Z-coordinate')
#GGraph.line( gp_adv, true, 'exchange'=>true, \
#'title'=>'CLOUD DENSITY TENDENCY', \
#'min'=> __split_range($OPT_range)[0], \
#'max'=> __split_range($OPT_range)[1], \
#'index'=>22 )  # red
#GGraph.line( gp_adv, false, 'exchange'=>true, 'index'=>22 ) # red
#GGraph.line( gp_diff, false, 'exchange'=>true, 'index'=>42 ) # blue
GGraph.line( gp_diff, false, 'exchange'=>true, 'index'=>102 ) # purple
GGraph.line( gp_others, false, 'exchange'=>true, 'index'=>42 ) # blue
#GGraph.line( gp_turb, false, 'exchange'=>true, 'index'=>92 ) # aqua
#GGraph.line( gp_fill, false, 'exchange'=>true, 'index'=>102 ) # purple
GGraph.line( gp_cond, false, 'exchange'=>true, 'index'=>32 ) # green
DCL.grcls