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

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

var = ARGV.shift
aomega = ARGV.shift

randoms = (101..125).to_a
nrandom = randoms.length

path = "../data/"
for n in 0...nrandom
  random = randoms[n]
p  fname = "#{path}Omega#{aomega}/baro_rn4expcs_freedecay_cpg100_Omega#{aomega}-#{random}.nc"

  file = NetCDF.open(fname)
  if var=="potvor"
    vor = file.var("vor")[true,true,-1]
  else
    ndat = file.var(var)
    ndat = ndat[true,true,-1]
  end

  if n==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]
    if var=="potvor"
      vname = "Potential Vorticity"
    else
      vname = file.var(var).att("long_name").get
    end
    dat = NArray.sfloat(lon.length,lat.length,nrandom)
  end
  file.close


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


  dat[true,true,n] = ndat
end


dat = dat.mean(2)

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)
    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
end

DCL::grcls