require "numru/ggraph"
include NumRu

#< interpret command-line arguments > 

#iws = (ARGV[0] || 1).to_i
iws = 4
#varname = ARGV[0].to_s

# variable settings

Grav = 8.87
Cp   = 820.0

varname = "Temp"
#varfile = varname+".nc"
varfile = 'venus_'+varname+".nc"
Hname = "Height"
#Hfile = Hname + ".nc"
Hfile = "venus_" + Hname + ".nc"

loncut=0
latcut=0
timecut=1000

#< open data > 

var = GPhys::IO.open(varfile,varname).cut('lon'=>loncut,'lat'=>latcut,'time'=>timecut)
height = GPhys::IO.open(Hfile,Hname).cut('lon'=>loncut,'lat'=>latcut,'time'=>timecut)

var.set_assoc_coords([height])
p "var:",var

#< prepare a height coordinate variable >

height_crd = VArray.new( 
   NArray[0.0, 0.5e4, 1e4, 1.5e4, 2e4, 2.5e4, 3e4, 3.5e4, 4e4, 4.5e4, 5e4, 5.5e4, 
          6e4, 6.5e4, 7e4, 7.5e4, 8e4, 8.5e4, 9e4, 9.5e4, 10.0e4],
   {"units"=>"m"}, Hname )

#< transform the vertical coordinate to theta >

var_on_height = var.interpolate("sig"=>height_crd)

km = height_crd.shape

temp=var_on_height.val
zh=height_crd.val

p temp, zh

vb = NArray.float(km[0]-1)

for k in 1..km[0]-1 do
  vb[k-1] = Grav/((temp[k]+temp[k-1])/2)* (Grav/Cp-(temp[k]-temp[k-1])/(zh[k]-zh[k-1]))
  print k, ' ', vb[k-1],"\n"
end



#< graphics >

#DCL.swpset('iwidth',800)
#DCL.swpset('iheight',400)
#DCL.swpset('ldump',true) if iws==4
#DCL.gropn(iws)
#DCL.sldiv('y',2,1)
#DCL.sgpset('isub', 96)      # control character of subscription: '_' --> '`'
#DCL.glpset('lmiss',true)

#DCL.grcls