# -*- coding: euc-jp -*-
# Title: Ruby script drawing contour map for deepconv/arare5 output data 
#
# History: 2011/09/27 (Masatsugu Odaka)
#
require "numru/ggraph"
include NumRu

filehead = "BS1998_"
filedir = "../time_000180000-000280000"
cuttime = 180000..266400
adv_file0   = "#{filedir}/#{filehead}PTempAdv_rank000000.nc"
adv_file1   = "#{filedir}/#{filehead}PTempAdv_rank000001.nc"
adv_file2   = "#{filedir}/#{filehead}PTempAdv_rank000002.nc"
adv_file3   = "#{filedir}/#{filehead}PTempAdv_rank000003.nc"
adv_file4   = "#{filedir}/#{filehead}PTempAdv_rank000004.nc"
adv_file5   = "#{filedir}/#{filehead}PTempAdv_rank000005.nc"
#nDiff = "#{filedir}/#{filehead}PTempNDiff.nc", 'PTempNDiff')
turb_file0  = "#{filedir}/#{filehead}PTempTurb_rank000000.nc"
turb_file1  = "#{filedir}/#{filehead}PTempTurb_rank000001.nc"
turb_file2  = "#{filedir}/#{filehead}PTempTurb_rank000002.nc"
turb_file3  = "#{filedir}/#{filehead}PTempTurb_rank000003.nc"
turb_file4  = "#{filedir}/#{filehead}PTempTurb_rank000004.nc"
turb_file5  = "#{filedir}/#{filehead}PTempTurb_rank000005.nc"
#disp  = "#{filedir}/#{filehead}PTempDisp.nc",  'PTempDisp' )
rad_file0   = "#{filedir}/#{filehead}PTempRad_rank000000.nc"
rad_file1   = "#{filedir}/#{filehead}PTempRad_rank000001.nc"
rad_file2   = "#{filedir}/#{filehead}PTempRad_rank000002.nc"
rad_file3   = "#{filedir}/#{filehead}PTempRad_rank000003.nc"
rad_file4   = "#{filedir}/#{filehead}PTempRad_rank000004.nc"
rad_file5   = "#{filedir}/#{filehead}PTempRad_rank000005.nc"

sfc_file0   = "#{filedir}/#{filehead}PTempSfc_rank000000.nc"
sfc_file1   = "#{filedir}/#{filehead}PTempSfc_rank000001.nc"
sfc_file2   = "#{filedir}/#{filehead}PTempSfc_rank000002.nc"
sfc_file3   = "#{filedir}/#{filehead}PTempSfc_rank000003.nc"
sfc_file4   = "#{filedir}/#{filehead}PTempSfc_rank000004.nc"
sfc_file5   = "#{filedir}/#{filehead}PTempSfc_rank000005.nc"

#adv = GPhys::IO.open(adv_file0, "PTempAdv" )
#turb = GPhys::IO.open(turb_file0, "PTempTurb" )
#rad = GPhys::IO.open(rad_file0, "PTempRad" )
#sfc = GPhys::IO.open(sfc_file0, "PTempSfc" )
#adv = GPhys::IO.open([adv_file0,adv_file1,adv_file2,adv_file3,adv_file4], "PTempAdv" )
#turb = GPhys::IO.open([turb_file0,turb_file1,turb_file2,turb_file3,turb_file4], "PTempTurb" )
#rad = GPhys::IO.open([rad_file0,rad_file1,rad_file2,rad_file3,rad_file4], "PTempRad" )
#sfc = GPhys::IO.open([sfc_file0,sfc_file1,sfc_file2,sfc_file3,sfc_file4], "PTempSfc" )
adv_a = GPhys::IO.open([adv_file0,adv_file1,adv_file2,adv_file3,adv_file4], "PTempAdv" )
turb_a = GPhys::IO.open([turb_file0,turb_file1,turb_file2,turb_file3,turb_file4], "PTempTurb" )
rad_a = GPhys::IO.open([rad_file0,rad_file1,rad_file2,rad_file3,rad_file4], "PTempRad" )
sfc_a = GPhys::IO.open([sfc_file0,sfc_file1,sfc_file2,sfc_file3,sfc_file4], "PTempSfc" )
dens0 = GPhys::IO.open("./restart_long-DensBZ.nc",   'DensBZ'  )
denscut = dens0
time  = GPhys::IO.open("#{filedir}/#{filehead}PTempAdv_rank000000.nc", 't').val
zz  = GPhys::IO.open("#{filedir}/#{filehead}PTempAdv_rank000000.nc", 'z').val

adv = adv_a.cut('t'=>cuttime)
turb = turb_a.cut('t'=>cuttime)
rad = rad_a.cut('t'=>cuttime)
sfc = sfc_a.cut('t'=>cuttime)

turb  = turb + sfc

# x-y mean

dz = 200.0
cpdry = 891.0
#=begin
adv_XYmean  = adv.average('x').average('y')
turb_XYmean = turb.average('x').average('y')
rad_XYmean  = rad.average('x').average('y')
#disp_XYmean  = disp.average('x').average('y')
#nDiff_XYmean  = nDiff.average('x').average('y')

# <Open DCL>
DCL.gropn(2)
# <Drawing data by GGraph>

#comax = 3.0e2
#comin = -3.0e2
comax = 3.0e-4
comin = -3.0e-4

total = adv_XYmean + turb_XYmean + rad_XYmean 
p total
GGraph.set_axes('xtitle'=>'Heating rate','ytitle'=>'z')
#total = total.mean('t')
#for i in 0...time.length do
  GGraph.line( total.mean('t'), \
#  GGraph.line( total.cut('t'=>time[i]), \
               true, 'exchange'=>true , \
               'index'=>12, 'type'=>1, 'label'=>'total', \
               'title'=>'Total PTemp Tendency ', \
               'max'=>comax, 'min'=>comin)    # °Üή
#end
#=begin
  GGraph.line( adv_XYmean.mean('t'), \
#  GGraph.line( adv_XYmean.cut('t'=>time[i]), \
               false, 'exchange'=>true , \
               'index'=>22, 'type'=>1, 'label'=>'adv', \
               'max'=>comax, 'min'=>comin)    # °Üή

  GGraph.line( turb_XYmean.mean('t'), \
#  GGraph.line( turb_XYmean.cut('t'=>time[i]), \
               false, 'exchange'=>true , \
                'index'=>32, 'type'=>1, 'label'=>'turb', \
               'max'=>comax, 'min'=>comin)    # Íðή³È»¶
  GGraph.line( rad_XYmean.mean('t'), \
#  GGraph.line( rad_XYmean.cut('t'=>time[i]), \
               false, 'exchange'=>true , \
               'index'=>42, 'type'=>1, 'label'=>'rad', \
               'max'=>comax, 'min'=>comin)    # Êü¼Í²ÃÇ®

#=end
=begin
  GGraph.line( disp_XYmean.mean('t'), \
               false, 'exchange'=>true , \
                'index'=>52, 'type'=>1, 'label'=>'disp', \
               'max'=>comax, 'min'=>comin)    # Íðή³È»¶

  GGraph.line( nDiff_XYmean.mean('t'), \
               false, 'exchange'=>true , \
                'index'=>62, 'type'=>1, 'label'=>'ndiff', \
               'max'=>comax, 'min'=>comin)    # Íðή³È»¶

=end
#  GGraph.line( - rad_int_XYmean.cut('t'=>time[i]), \
#end
# <Close DCL>
DCL.grcls