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

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


def max(val)
  tmp = (Math::log(val)/Math::log(10)).floor
  tmp = 10.0**tmp
  val /= tmp
  if val>7
    tmp *= 10
  elsif val>3
    tmp *= 5
  elsif val>2.1
    tmp *= 2.5
  elsif val>1.2
    tmp *= 2
  end
  return tmp
end

fname = ARGV.shift
file = NetCDF.open(fname)
vor = file.var("vor")[true,true,0]
strfunc = file.var("strfunc")[true,true,0]
lon = file.var("lon").get
lat = file.var("lat").get
t = file.var("t")[-1][0]
omega = file.att("Omega").get[0]
a = file.att("Radius").get[0]
im,jm = vor.shape
nm = file.dim("n").length
file.close

w_Initial(nm,im,jm)

f = 2*omega*NMath::sin(lat*Math::PI/180)
#povor = (vor+f.reshape!(1,jm))/(hsfc+h)

energy = n_EnergyFromStreamFunc_w(w_xy(strfunc))
enstrophy = n_EnstrophyFromStreamFunc_w(w_xy(strfunc))

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",14)
DCL::gropn(iws)

DCL::uzfact(0.8)
DCL::uwsgxa(lon)
DCL::uwsgya(lat)

#dat = [vor,povor]
dat = [vor]
title = ["Vorticity","Potential Vorticity"]
for p in 0...dat.length
  DCL::grfrm
  DCL::grsvpt(0.1,0.95,0.08,0.65)
  DCL::grswnd(0,360,-90,90)
  DCL::grstrn(1)
  DCL::grstrf

  ntone = 10
  tmp = dat[p].abs.max
  max = max(tmp)
  ci = max/10
  DCL::ueitlv
  for i in 0...ntone
    DCL::uestlv(-max*(i+1)/ntone,-max*i/ntone,(57-5*i)*1000+999)
    DCL::uestlv(max*(i)/ntone,max*(i+1)/ntone,(52+5*i)*1000+999)
  end
  DCL::uestlv(-999,-max,10999)
  DCL::uestlv(max,999,99999)


  DCL::uetone(dat[p])
  if ci >= 0.5
    DCL::udsfmt("(i8)")
  end
  DCL::udgcla(-ci*20,ci*20,ci)
  DCL::udcntz(dat[p])

  DCL::usdaxs
  DCL::uxsttl("t","#{title[p]} t=0",0)
  DCL::uxsttl("b","longitude",0)
  DCL::uysttl("l","latitude",0)

  DCL::sgtxzr(0.95,0.015,"Contour interval=#{ci}",0.012,0,1,1)
end

dat = [energy,enstrophy]
title = ["Energy", "Enstrophy"]
DCL::sglset("lclip",true)
for i in 0...2
  DCL::grfrm
  DCL::grsvpt(0.1,0.95,0.08,0.65)
  DCL::grswnd(0,nm,[1.0e-10,1.0e-5][i],[1.0e-1,1.0e2][i])
  DCL::grstrn(2)
  DCL::grstrf
  mask = dat[i].gt(0)
  DCL::sgplzu(NArray.sfloat(nm+1).indgen[mask],dat[i][mask],1,2)
  DCL::usdaxs
  DCL::uxsttl("t","#{title[i]} spectrum t=0",0)
  DCL::uxsttl("b","wave number",0)
end

DCL::grcls