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

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

#ntim = 0
ntim = -1

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)


iws = 4
if ARGV.delete("-ps")
  iws = 2
elsif ARGV.delete("-png")
  DCL::swlset("lwnd",false)
end

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

DCL::uzfact(0.7)

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



DCL::uzrset("uxuser",0)
  u = var_u[true,true,ntim].val
  hsfc = var_hsfc[true,true,ntim].val if var_hsfc
  vor = var_vor[true,true,ntim].val
  ux = u.mul_add(lonw,0) 
  ucos = ux*cosphi
  if var_hsfc
#    uhcos = (u*(hsfc+h)).mul_add(lonw,0)*cosphi
    gh = (g*hsfc).mul_add(lonw,0)
    tmp = (vor+f)/(hsfc+h)
  else
#    uhcos = nil
    gh = 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]
  data = [ux, ucos, gh, 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
    when 1
      xmax = 0.5
    when 2
      if data[n]
        xmax = data[n].abs.max
        if xmax == 0
          xmax = 1
        else
          xmax = Math::log(xmax)/Math::log(10)
          xmax = 10**(xmax.ceil)
          if xmax/data[n].abs.max>5
            xmax = xmax/5
          elsif xmax/data[n].abs.max>2
            xmax = xmax/2
          end
        end
      else
        xmax = 1
      end
    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[ntim]}sec",0.02,0,-1,2)
  if t[ntim]>0
    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