require "numru/gphys"
require "spml"
require "numru/dcl"
include NumRu
include SPML

ARGV.length<1 && raise("Usage: ruby #$0 fname")

rmiss = -9999.9



fname = ARGV.shift
vname = "vor"
file = NetCDF.open(fname)
dat_all = file.var(vname)
lon = file.var("lon").get
lat = file.var("lat").get
t = file.var("t").get
omega = file.att("Omega").get[0]
g = file.att("Grav")
g = g.get[0] if g
h = file.att("Hbar")
h = h.get[0] if h
a = file.att("Radius").get[0]
im = lon.length
jm = lat.length
ntim = t.length
nm = file.dim("n").length

w_Initial(nm,im,jm)

f = 2*omega*NMath::sin(lat*Math::PI/180)

iws = ARGV.delete("-ps") ? 2 : 4
if ARGV.delete("-png")
  DCL::swlset("lwait",false)
  DCL::swlset("ldump",true)
  DCL::swlset("lwnd",false)
end
DCL::sglset("lfull",true)
DCL::udlset("lmsg",false)
DCL::swiset("iclrmap",5)
#DCL::swiset("iclrmap",14)
DCL::gllset("lmiss",true)
DCL::glrset("rmiss",rmiss)
DCL::uzlset("labelxt",true)
DCL::gropn(iws)


DCL::uzfact(0.8)

  ntone = 6
  fact = 0.25
  tone =( NArray.sfloat(ntone+1).indgen-ntone )*fact
  for i in 0...ntone
    DCL::uestlv(tone[i],tone[i+1],(99-i*(99/ntone).to_i)*1000+999)
#    DCL::uestlv(tone[i],tone[i+1],(55+i*(44.0/ntone).to_i)*1000+999)
  end
  DCL::uestlv(tone.max,999,99999)

for nt in [0]
#for nt in [-1]
#for nt in 0...ntim
spect = w_xy(dat_all[true,true,nt])
nspect = spect.length
dat = NArray.sfloat(nm+1,nm+1).fill(rmiss)

for l in 1..nspect
  n,m = ISPACK.snl2nm(l)
  m = m.abs
  if dat[m-1,n-1]>0
    dat[m-1,n-1] = Math::sqrt( dat[m-1,n-1]**2 + spect[l-1]**2 )
  else
    dat[m-1,n-1] = spect[l-1].abs
  end
end

log10 = Math::log(10)
for n in 0..nm
  for m in 0..nm
    if dat[m,n]>0
      dat[m,n] = Math::log(dat[m,n])/log10
    end
  end
end


  DCL::grfrm
  DCL::grsvpt(0.2,0.8,0.02,0.62)
  DCL::grswnd(0,30,0,30)
#  DCL::grswnd(0,100,0,100)
#  DCL::grswnd(0,nm,0,nm)
  DCL::grstrn(1)
  DCL::grstrf

DCL::uwsgxa(NArray.sint(nm+1).indgen)
DCL::uwsgya(NArray.sint(nm+1).indgen)

  DCL::uetonf(dat)

#  DCL::usdaxs
  DCL::uxaxdv("t",5,10)
  DCL::uyaxdv("l",5,10)
  DCL::uysttl("l","total wave number",0)
  DCL::uxsttl("t","zonal wave number",0)

  DCL::sgtxzr(0.75,0.15,"t=#{"%.1f"%t[nt]}",0.02,0,1,1)


  DCL::grfig
  DCL::grswnd(tone.min,tone.max,0,1)
  DCL::grsvpt(0.5,0.75,0.1,0.12)
  DCL::grstrn(1)
  DCL::grstrf
DCL::uwsgxa(tone)
DCL::uwsgya([0,1])

  DCL::uetone(NArray[tone,tone])
  DCL::uzfact(0.5)
  DCL::uxaxdv("b",fact,0.5)
  DCL::uxsttl("b","log_10\"",0)
  DCL::uzfact(2)

end


DCL::grcls

file.close