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

filebase = 'restart_long-PTempBZ.nc'
filedir = "time_000180000-000280000"
filehead = "BS1998_PTemp"
file0 = "../#{filedir}/#{filehead}_rank000000.nc"
file1 = "../#{filedir}/#{filehead}_rank000001.nc"
file2 = "../#{filedir}/#{filehead}_rank000002.nc"
file3 = "../#{filedir}/#{filehead}_rank000003.nc"
file4 = "../#{filedir}/#{filehead}_rank000004.nc"
file5 = "../#{filedir}/#{filehead}_rank000005.nc"
varname1 = 'PTemp'
var1  = GPhys::IO.open([file0,file1,file2,file3,file4,file5], varname1)
time  = GPhys::IO.open(file1, 't').val

varname0 = 'PTempBZ'

#comax = 249.0
#comin = 247.0
comax = 0.010
comin = -0.0010
var0  = GPhys::IO.open(filebase, varname0)
dptdz = GPhys::IO.open('../DPTDz.nc', 'DPTDz')
dptdz = dptdz.cut('z'=>300..20000)

time  = ARGV[0].to_f
timelabel = time.to_s
x_co  = GPhys::IO.open(file1, 'x').val

# <Open DCL>
DCL.swlset( 'ldump', true )
DCL.gropn(2)

varall = var1.cut('t'=>180000..266400).mean('t').cut(true,true,true) + var0

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

varmean = varall.average('x').average('y')
stability = varmean.deriv('z')

GGraph.line( dptdz, 
             true, 'exchange'=>true ,
             'index'=>2, 'type'=>1, 'label'=>'init',
             'title'=>'Stability.',
             'max'=>comax, 'min'=>comin)    
#=begin
GGraph.line( stability, 
             false, 'exchange'=>true ,
             'index'=>1, 'type'=>3, #'label'=>'after',
             'max'=>comax, 'min'=>comin)    
#=end

# <Close DCL>
DCL.grcls