# -*- 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-all_"

adv   = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempAdv.nc",   'PTempAdv'  )
nDiff = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempNDiff.nc", 'PTempNDiff')
turb  = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempTurb.nc",  'PTempTurb' )
disp  = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempDisp.nc",  'PTempDisp' )
rad   = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempRad.nc",   'PTempRad'  )
sfc   = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempSfc.nc",   'PTempSfc'  )
dens0 = GPhys::IO.open("./restart_long-DensBZ.nc",   'DensBZ'  )
exner0 = GPhys::IO.open("./restart_long-ExnerBZ.nc",   'ExnerBZ'  )

#time  = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempAdv.nc", 't').val

zz  = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempAdv.nc", 'z').val

turb  = turb + sfc
cuttime = 80000
nz      = 100
dz      = 200.0
cpdry   = 891.0
zlength = 20000.0
radinit = - 293.1
turbinit = 189.70
#turbinit = 293.1
comax   = 500.0
comin   = -500.0

#x-y 平均し, K.s-1 を W.m-3 に変換する
denscut = dens0.mean('x').mean('y')
adv_XYmean =  adv.mean('x').mean('y').cut('t'=>cuttime) * cpdry * denscut #* exner0.mean('x').mean('y')
adv_XYmean.units = 'W.m-2'
turb_XYmean =  turb.mean('x').mean('y').cut('t'=>cuttime) * cpdry * denscut* exner0.mean('y').mean('x')
#turb_XYmean = - turb.mean('x').mean('y').mean('t') * cpdry * denscut * exner0.mean('y').mean('x')
#turb = - turb.cut('t'=>cuttime) * cpdry * dens0
#turb_XYmean = turb.mean('x').mean('y')
#turb_XYmean = - turb.mean('x').mean('y').cut('t'=>cuttime) * cpdry * denscut
turb_XYmean.units = 'W.m-2'
rad_XYmean  = rad.mean('x').mean('y').cut('t'=>cuttime) * cpdry * denscut * exner0.mean('y').mean('x')
rad_XYmean.units = 'W.m-2'
disp_XYmean  = disp.mean('x').mean('y').cut('t'=>cuttime) * cpdry * denscut
disp_XYmean.units = 'W.m-2'
nDiff_XYmean  = nDiff.mean('x').mean('y').cut('t'=>cuttime) * cpdry * denscut
nDiff_XYmean.units = 'W.m-2'


# 鉛直方向に積分し, W.m-2 を求める
#turb_XYmean[zz.length - 1] = turbinit
turb_XYmean[0] = turbinit
rad_XYmean[zz.length - 1] = radinit
for i in 1...zz.length do
  adv_XYmean[i] = adv_XYmean[i-1].val - adv_XYmean[i].val * dz 
  turb_XYmean[i] = turb_XYmean[i-1].val - turb_XYmean[i].val * dz 
  disp_XYmean[i] = disp_XYmean[i-1].val - disp_XYmean[i].val * dz 
  nDiff_XYmean[i] = nDiff_XYmean[i-1].val - nDiff_XYmean[i].val * dz 
  # 太陽放射フラックスは上端から積分
  j = zz.length - i 
  rad_XYmean[j-1] = rad_XYmean[j].val + rad_XYmean[j-1].val * dz 
#  turb_XYmean[j-1] = turb_XYmean[j].val + turb_XYmean[j-1].val * dz 
end
DCL.gropn(2)

# <Drawing data by GGraph>

  GGraph.line( adv_XYmean, \
               true, 'exchange'=>true , \
               'index'=>22, 'type'=>1, 'label'=>'adv', \
               'title'=>'Energy Flux', \
               'max'=>comax, 'min'=>comin)    # 移流

  GGraph.line( turb_XYmean, \
               false, 'exchange'=>true , \
                'index'=>32, 'type'=>1, 'label'=>'turb', \
               'max'=>comax, 'min'=>comin)    # 乱流拡散

  GGraph.line( disp_XYmean, \
               false, 'exchange'=>true , \
                'index'=>52, 'type'=>1, 'label'=>'disp', \
               'max'=>comax, 'min'=>comin)    # 乱流拡散

  GGraph.line( nDiff_XYmean, \
               false, 'exchange'=>true , \
                'index'=>62, 'type'=>1, 'label'=>'ndiff', \
               'max'=>comax, 'min'=>comin)    # 乱流拡散

  GGraph.line( rad_XYmean, \
               false, 'exchange'=>true , \
               'index'=>42, 'type'=>1, 'label'=>'rad', \
               'max'=>comax, 'min'=>comin)    # 放射加熱

# <Close DCL>
DCL.grcls