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

ARGV.length<4 && raise("Usage: ruby #$0 var Omega G U [t]")

var = ARGV.shift
omega = ARGV.shift
g = ARGV.shift
u = ARGV.shift
t0 = ARGV[0]
if t0.to_i.to_s == t0
  ARGV.shift
  t0 = t0.to_i
else
  t0 = nil
end


fname = "../data/spshallow-zd_rn4cn_freedecay_Omega#{omega}_G#{g}_U#{u}-2.nc"
if !File.exist?(fname)
  fname = "../data/spshallow-zd_rn4cn_freedecay_Omega#{omega}_G#{g}_U#{u}.nc"
end

p fname
file = NetCDF.open(fname)
t = file.var("t").get
if t0
  tmp = (t-t0).abs
  ntim = tmp.eq(t.min).where[0]
else
  ntim = -1
end
t = t[ntim]

if var=="potvor"
  vor = file.var("vor")[true,true,ntime]
  hsfc = file.var("hsfc")[true,true,ntim]
  vname = "Potential Vorticity"
else
  dat = file.var(var)
  vname = dat.att("long_name").get
  dat = dat[true,true,ntim]
end
lon = file.var("lon").get
lat = file.var("lat").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]
file.close


if var=="potvor"
  f = 2*omega*NMath::sin(lat*Math::PI/180)
  dat = (vor+f.reshape(1,lat.length))/(hsfc+h)
end

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::gropn(iws)

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




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

ci = max/10
if var=="vor"
  if t==0
    ci = 100
  elsif t==10
    ci = 20
  elsif t==5
    if (omega==0||omega==25||omega==100) && g==10
      ci = 50
    end
  end
elsif var=="strfunc"
  if t==0
    ci = 0.05
  elsif t==5
    if omega==25 && g==10
      ci = 0.1
    end
  end
end


if ci >= 0.5
  DCL::udsfmt("(i8)")
end
DCL::udgcla(-ci*20,ci*20,ci)
DCL::uddclv(0)
=begin
unless var=="div"
  DCL::udsclv(0,2,1,"",0)
end
=end

=begin
max *= 0.5
ntone = 5
DCL::uestlv(-max/ntone,max/ntone,60999)
for i in 1...ntone
  DCL::uestlv(-max*(i+1)/ntone,-max*i/ntone,(60-(36.0/(ntone-1)).to_i*i)*1000+999)
  DCL::uestlv(max*(i)/ntone,max*(i+1)/ntone,(60+(36.0/(ntone-1)).to_i*i)*1000+999)
end
DCL::uestlv(-999e5,-max,20999)
DCL::uestlv(max,999e5,99999)
=end
DCL::uestlv(-999e5,-max/2,45)
DCL::uestlv(-max/2,0,34)


#DCL::umlset("lgridmn",false)


      DCL::grfrm
      DCL::grsvpt(0.0,1.0,0.005,0.715)
    DCL::grsmpl(90,0,0)
    DCL::grstrn(30)
    DCL::umpfit
    DCL::grstrf

      DCL::uetone(dat)
      #DCL::uetonf(dat)
    DCL::udcntz(dat)

    DCL::umpglb

    DCL::sgtxzr(0.02,0.68,"#{vname}",0.02,0,-1,2)
    DCL::sgtxzr(0.02,0.65,"t=#{t}sec",0.02,0,-1,2)
    unless t==0
      if g
        DCL::sgtxzr(0.98,0.035,"Omega=#{omega.to_i}, g=#{g.to_i}, U=#{u}, a=#{a.to_i}",0.012,0,1,1)
      else
        DCL::sgtxzr(0.98,0.035,"Omega=#{omega.to_i}, g=inf, U=#{u}, a=#{a.to_i}",0.012,0,1,1)
      end
    end
    DCL::sgtxzr(0.98,0.015,"Contour interval=#{ci}",0.012,0,1,1)

DCL::grcls