require "numru/ggraph"
include NumRu

ZMAX=1.0e04
ZMIN=0.0
VERTICAL_LVNUM=20
OUTPUT_TICK=7200.0
IC_RADIUS=6371e03

fieldz_fileName="intr_fieldzData.nc"
metric_fileName="intr_TFGridData.nc"
cuts_hash = { 'lon'=>0 }

xyz_Phys = GPhys::NetCDF_IO.open(fieldz_fileName, 'xyz_R')
xy_Zs = GPhys::NetCDF_IO.open(metric_fileName, 'xy_Zs').cut('time'=>0)
xy_G_1div2 = GPhys::NetCDF_IO.open(metric_fileName, 'G_1div2').cut('time'=>0)

tlevelNum = xyz_Phys.shape[3]
levelNum = xyz_Phys.shape[2]
latNum = xyz_Phys.shape[1]
lonNum = xyz_Phys.shape[0]

Xiz = xyz_Phys.axis("levels").to_gphys
grid = xyz_Phys.cut("time"=>0).grid_copy
xyz_Z = GPhys.new(grid, VArray.new(NArray.float(lonNum,latNum,levelNum), {"units"=>"m"}, "xyz_Z"))
xyz_Zs = GPhys.new(grid, VArray.new(NArray.float(lonNum,latNum,levelNum), {"units"=>"m"}, "xyz_Zs"))
xyz_G_1div2 = GPhys.new(grid, VArray.new(NArray.float(lonNum,latNum,levelNum), {"units"=>"m"}, "xy_G_1div2"))

for levelID in 0..levelNum-1
  xyz_Zs.data[true, true, levelID] = xy_Zs.data[true, true]
  xyz_G_1div2.data[true, true, levelID] = xy_G_1div2.data[true, true]
end

xyz_Z = Xiz * ( 1.0 - xyz_Zs / ZMAX ) + xyz_Zs 
xyz_Z.name = "z"
xyz_Z.long_name = "altitude"

xyz_Phys.set_assoc_coords([xyz_Z])

p "xyz_Phys:", xyz_Phys

delZ = ( ZMAX - ZMIN ) / VERTICAL_LVNUM 
z = NArray.float(VERTICAL_LVNUM).indgen(ZMIN, delZ)

z_crd = VArray.new(z, {"units"=>"m"}, "z")


DCL.gropn(1)
DCL.sgpset('lcntl', false)
DCL.uzfact(0.7)

for tlevel in 0..tlevelNum-1

  time=OUTPUT_TICK*tlevel

  xyz_Phys_onZ_slice = xyz_Phys.interpolate("levels"=> z_crd).cut(cuts_hash).cut("time"=>time)
  xy_G_1div2_slice = xyz_G_1div2.cut(cuts_hash)
  gamma_slice = 1.0 + xyz_Z.cut(cuts_hash) / IC_RADIUS

  xyz_Var_slice = xyz_Phys_onZ_slice / ( xy_G_1div2_slice * gamma_slice**2 )   

  GGraph.set_fig('viewport'=>[0.1,0.9,0.3,0.8])
  GGraph.tone(xyz_Var_slice)
  GGraph.contour(xyz_Var_slice, false)
  GGraph.color_bar('vlength'=>0.5,"landscape"=>true,'tickintv'=>0)

end

DCL.grcls