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

ARGV.length==5 || raise("Usage: ruby #$0 var Omega G U map")

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

path = "../data/"
fname1 = "#{path}spshallow-zd_rn4cn_freedecay_Omega#{omega}_G#{g}_U#{u}.nc"
fname2 = "#{path}spshallow-zd_rn4cn_freedecay_Omega#{omega}_G#{g}_U#{u}-2.nc"

file1 = NetCDF.open(fname1)
file = [file1]
if File.exist?(fname2)
  file2 = NetCDF.open(fname2)
  file = [file1,file2]
else
  file = file1
end

if var=="potvor"
  vor = GPhys::IO.open(file,"vor")
  hsfc = GPhys::IO.open(file,"hsfc")
  vname = "Potential Vorticity"
  lon = vor.coord("lon").val
  lat = vor.coord("lat").val
  t = vor.coord("t").val
else
  dat = GPhys::IO.open(file,var)
  vname = dat.get_att("long_name")
  lon = dat.coord("lon").val
  lat = dat.coord("lat").val
  t = dat.coord("t").val
end
omega = file1.att("Omega").get[0]
g = file1.att("Grav")
g = g.get[0] if g
h = file1.att("Hbar")
h = h.get[0] if h
a = file1.att("Radius").get[0]



DCL::swlset("lwait",false)
DCL::swlset("ldump",true)
DCL::swlset("lwnd",false)
DCL::sglset("lfull",true)
DCL::udlset("lmsg",false)
DCL::swiset("iclrmap",14)
DCL::swiset("iwidth",450)
DCL::swiset("iheight",325)
DCL::gropn(4)

DCL::uzfact(0.7)



if var=="potvor"
  f = 2*omega*NMath::sin(lat*Math::PI/180).reshape!(1,lat.length)
  tmp = (vor[true,true,-1].val+f)/(hsfc[true,true,-1].val+h)
  tmp = tmp.abs.max
else
  tmp = dat[true,true,-1].val.abs.max
end

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)
tone = NArray.sfloat(ntone*2+1,2)
tone[true,0] = (NArray.sfloat(ntone*2+1).indgen/ntone-1)*max
tone[true,1] = tone[true,0]

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)

DCL::uzlset("labelyl",false)


if map=="c"
  nfig = 1
  itr = 10
elsif map=="o" || map=="p"
  nfig = 2
  itr = 30
end

t.length.times do |m|
  if var=="potvor"
    tmp = (vor[true,true,m].val+f)/(hsfc[true,true,m].val+h)
  else
    tmp = dat[true,true,m].val
  end
  DCL::uwsgxa(lon)
  DCL::uwsgya(lat)
  nfig.times do |n|
    if n==0
      DCL::grfrm
    else
      DCL::grfig
    end
    if map=="c"
      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 map=="c"
      DCL::grsmpl(180,90,0)
    elsif map=="o"
      if n==0
        DCL::grsmpl(90,0,0)
      else
        DCL::grsmpl(270,0,0)
      end
    elsif map=="p"
      if n==0
        DCL::grsmpl(90,90,0)
      else
        DCL::grsmpl(270,-90,0)
      end
    end
    DCL::grstrn(itr)
    DCL::umpfit
    DCL::grstrf

    if map=="p"
      DCL::uetone(tmp)
    else
      DCL::uetonf(tmp)
    end
    #DCL::udcntz(tmp)

    DCL::umpglb

    DCL::sgtxzr(0.02,0.68,"#{vname} t=#{"%.1f"%t[m]}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

  DCL::uwsgxa(tone[true,0])
  DCL::uwsgya([0,1])
  DCL::grfig
  DCL::grsvpt(0.75,0.95,0.08,0.1)
  DCL::grswnd(-max,max,0,1)
  DCL::grstrn(1)
  DCL::grstrf
  DCL::uetone(tone)
  DCL::usdaxs
end

DCL::grcls