require "gtk2"
require "numru/gphys"
require "numru/dcl"
require "spml"
include NumRu
include SPML

ARGV.length==3 || raise("Usage: ruby #$0 omega g u0")

omega = ARGV.shift
g = ARGV.shift
u0 = ARGV.shift

path = "../data/"
fname1 = "#{path}spshallow-zd_rn4cn_freedecay_Omega#{omega}_G#{g}_U#{u0}.nc"
fname2 = "#{path}spshallow-zd_rn4cn_freedecay_Omega#{omega}_G#{g}_U#{u0}-2.nc"
if File.exist?(fname2)
  file1 = NetCDF.open(fname1)
  file2 = NetCDF.open(fname2)
  file = [file1,file2]
else
  file1 = NetCDF.open(fname1) 
  file = file1
end

var_u = GPhys::IO.open(file,"vellon")
var_vor = GPhys::IO.open(file,"vor")
var_hsfc = file1.var("hsfc") ? GPhys::IO.open(file,"hsfc") : nil
t = GPhys::IO.open(file,"t").val
lon = file1.var("lon").get
lat = file1.var("lat").get
lonw = file1.var("lon_weight")
if lonw
  lonw = lonw.get
else
  lonw = NArray.sfloat(lon.length).fill(1)
end
lonw /= lonw.sum
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]
nm = file1.dim("n").length

w_Initial(nm,lon.length,lat.length)

f = 2*omega*NMath::sin(lat*Math::PI/180)
f.reshape!(1,lat.length)
cosphi = NMath::cos(lat*Math::PI/180)

ntim = t.length

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

DCL::uzfact(0.7)

title = ["u", "ucos#{DCL::csgi(172)}", "uhcos#{DCL::csgi(172)}", "PV", 'PV_y"']



DCL::uzrset("uxuser",0)
for i in 0...ntim
#for i in [ntim-1]
  u = var_u[true,true,i].val
  hsfc = var_hsfc[true,true,i].val if var_hsfc
  vor = var_vor[true,true,i].val
  ux = u.mul_add(lonw,0) 
  ucos = ux*cosphi
  if var_hsfc
    uhcos = (u*(hsfc+h)).mul_add(lonw,0)*cosphi
    tmp = (vor+f)/(hsfc+h)
  else
    uhcos = nil
    tmp = vor+f
  end
  potvor = tmp.mul_add(lonw,0)
  potvor_y = (xy_GradLat_w(w_xy(tmp))/a).mul_add(lonw,0)

  data = [ux, ucos, uhcos, potvor, potvor_y]

  vymin = 0.13
  vymax = 0.66
  for n in 0...5
    n==0 ? DCL::grfrm : DCL::grfig
    vxmin = 0.02+0.195*n
    vxmax = vxmin+0.16
    case n
    when 0
      xmax = 1
    when 1,2
      xmax = 0.5
    when 3
      if omega==0
        xmax = 20
      else
        xmax = 2*omega
      end
    when 4
      if omega==0
        xmax = 1000
      else
        xmax = 2*omega
      end
    end
    DCL::grsvpt(vxmin,vxmax,vymin,vymax)
    DCL::grswnd(-xmax,xmax,-90,90)
    DCL::grstrn(1)
    DCL::grstrf

    DCL::usxaxs("b")
    DCL::uyaxdv("u",10,30)
    DCL::uxsttl("b",title[n],0)

    DCL::sgplzu(data[n],lat,1,2) if data[n]
  end
  DCL::sgtxzr(0.02,0.71,"t=#{"%3.1f"%t[i]}sec",0.02,0,-1,2)
  if g
    DCL::sgtxzr(0.98,0.007,"Omega=#{omega.to_i}, g=#{g.to_i}, H=#{h.to_i}, U=#{u0}, a=#{a.to_i}",0.012,0,1,1)
  else
    DCL::sgtxzr(0.98,0.007,"Omega=#{omega.to_i}, g=inf, U=#{u0}, a=#{a.to_i}",0.012,0,1,1)
  end

end

DCL::grcls