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

module NumRu
  class GPhys
    def deriv(dim)
      dim = dim_index(dim)    # input dim can be a String or an Integer
      x = axis(dim).to_gphys
      a = [true]*dim + [1..-1,false]
      b = [true]*dim + [0..-2,false]
      dydx = ( self[*a] - self[*b] ) / ( x[1..-1] - x[0..-2] )
      xi = ( x[1..-1] + x[0..-2] ) / 2
      dydx.axis(dim).set_pos(xi.coord(0))
      dydx.long_name = "d #{name}/d #{x.name}"
      dydx
    end
  end
end


include NumRu
filehead = "BS1998-all_"
file_rho = './restart_long-DensBZ.nc'
file1 = '../time_000000000-000172000/BS1998-all_PTemp.nc'
turbs  = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempTurb.nc",  'PTempTurb' )
sfc   = GPhys::IO.open("../time_000000000-000172000/#{filehead}PTempSfc.nc",   'PTempSfc'  )
varname_rho = 'DensBZ'
varname1 = 'PTemp'


#comax = 249.0
#comin = 247.0
comax = 400.0
comin = -400.0
#comax = 5.0e-4
#comin= -5.0e-4
kappa = 155.0
cuttime = 80000
var_rho  = GPhys::IO.open(file_rho, varname_rho)
var1  = GPhys::IO.open(file1, varname1)
x_co  = GPhys::IO.open(file1, 'x').val

turbs  = turbs + sfc
turb_XYmean = turbs.average('x').average('y').mean('t')


# <Open DCL>
DCL.gropn(2)

#varall = var1.mean('t').cut(true,true,true) * kappa
varall = var1.mean('x').mean('y').mean('t')


# 安定度の計算 dptdz
#stability = varall.cut(1,true,true,time)

var_rho_mean = var_rho.mean('x').mean('y')

dptdz = varall.deriv('z') * kappa * var_rho_mean[1..99] * 891.0
#dptdz = dptdz.deriv('z')
p dptdz

#varks = - dptdz.deriv('z')
#p turb
#for i in 1..97
#  p turb[i].val
#end
#p varks
p turb_XYmean
GGraph.line( - dptdz, 
#GGraph.line( varmean, 
             true, 'exchange'=>true ,
             'index'=>2, 'type'=>1, 'label'=>'init',
             'title'=>'Mean Potemtial Temp.',
             'max'=>comax, 'min'=>comin)    
GGraph.line( turb_XYmean[7..-1], 
#GGraph.line( varmean, 
             false, 'exchange'=>true ,
             'index'=>22, 'type'=>1, 'label'=>'turb',
             'max'=>comax, 'min'=>comin)    


# <Close DCL>
DCL.grcls