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

ARGV.length==4 || raise("Usage: ruby #$0 var Omega Random map")

var = ARGV.shift
omega = ARGV.shift
random = ARGV.shift
map = ARGV.shift

path = "../data/"
fname = "#{path}Omega#{omega}/baro_rn4expcs_freedecay_cpg100_Omega#{omega}-#{random}.nc"

file = NetCDF.open(fname)

if var=="potvor"
  vor = GPhys::IO.open(file,"vor")
  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 = file.att("Omega").get[0]
a = file.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
  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
  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)
    DCL::sgtxzr(0.98,0.035,"Omega=#{omega.to_i}, U=sqrt2, a=#{a.to_i}, random=#{random}",0.012,0,1,1)
    #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