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

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

var = ARGV.shift
omega = ARGV.shift
g = ARGV.shift
u = ARGV.shift

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

file = NetCDF.open(fname)
if var=="potvor"
  vor = file.var("vor")[true,true,-1]
  hsfc = file.var("hsfc")[true,true,-1]
  vname = "Potential Vorticity"
else
  dat = file.var(var)
  vname = dat.att("long_name").get
  dat = dat[true,true,-1]
end
lon = file.var("lon").get
lat = file.var("lat").get
t = file.var("t")[-1][0]
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",14)
DCL::gropn(iws)

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




p 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
ntone = 10
ci = max/10
for i in 1...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(-999e5,-max,10999)
DCL::uestlv(max,999e5,99999)

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

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


for nfig in 0...3
  [1,2,2][nfig].times do |n|
    if n==0
      DCL::grfrm
    else
      DCL::grfig
    end
    if nfig==0
      DCL::grsvpt(0.005,0.995,0.005,0.715)
    elsif n==0
      DCL::grsvpt(0.005,0.5,0.005,0.715)
    elsif n==1
      DCL::grsvpt(0.5,0.995,0.005,0.715)
    end
    if nfig==0
      DCL::grsmpl(180,90,0)
    elsif nfig==1
      if n==0
        DCL::grsmpl(90,0,0)
      else
        DCL::grsmpl(270,0,0)
      end
    elsif nfig==2
      if n==0
        DCL::grsmpl(90,90,0)
      else
        DCL::grsmpl(270,-90,0)
      end
    end
    DCL::grstrn([10,30,30][nfig])
    DCL::umpfit
    DCL::grstrf

    if nfig==2
      DCL::uetone(dat)
    else
      DCL::uetonf(dat)
    end
    DCL::udcntz(dat)

    DCL::umpglb

    DCL::sgtxzr(0.02,0.68,"#{vname} t=#{t}sec",0.02,0,-1,2)
    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
    DCL::sgtxzr(0.98,0.015,"Contour interval=#{ci}",0.012,0,1,1)
  end
end

DCL::grcls