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

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

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

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


u = GPhys::IO.open(file,"vellon")[true,true,1..-1]
hsfc = GPhys::IO.open(file,"hsfc")[true,true,1..-1] if var=~/h/
nlon,nlat,ntim = u.shape
lat = u.coord("lat").val
tim = u.coord("t").val
nm = file1.dim("n").length
lonw = file1.var("lon_weight")
lonw = lonw ? lonw.get : NArray.sfloat(nlon).fill(1)
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]


lonw /= lonw.sum
ux = NArray.sfloat(ntim,nlat)
coslat = NMath::cos(lat*Math::PI/180)
case var
when "u-xm"
  for n in 0...ntim
    ux[n,true] = u[true,true,n].val.mul_add(lonw,0)
  end
  umax = 2.0
when "ucos-xm"
  for n in 0...ntim
    ux[n,true] = u[true,true,n].val.mul_add(lonw,0)*coslat
  end
  umax = 0.5
when "uhcos-xm"
  for n in 0...ntim
    ux[n,true] = (u[true,true,n].val*(hsfc[true,true,n].val+h)).mul_add(lonw,0)*coslat
  end
  umax = 0.5
end


iws = ARGV.delete("-ps") ? 2 : 4
if ARGV.delete("-png")
  DCL::swlset("lwnd",false)
  DCL::swlset("ldump",true)
end
DCL::sglset("lfull",true)
DCL::udlset("lmsg",false)
DCL::swiset("iclrmap",14)
DCL::gropn(iws)

DCL::grfrm
DCL::grsvpt(0.09,0.98,0.08,0.67)
DCL::grswnd(0.1,tim[-1],-90,90)
DCL::grstrn(3)
DCL::grstrf
DCL::uzfact(0.7)

DCL::uwsgxa(tim)
DCL::uwsgya(lat)

ntone = 10
=begin
for i in 1...ntone
  DCL::uestlv(-umax*(i+1)/ntone,-umax*i/ntone,(57-5*i)*1000+999)
  DCL::uestlv(umax*(i)/ntone,umax*(i+1)/ntone,(52+5*i)*1000+999)
end
DCL::uestlv(-999,-umax,10999)
DCL::uestlv(umax,999,99999)
DCL::uetonf(ux)
=end
#DCL::uestlv(-999,-0.1,15)
DCL::uestlv(-999,0,15)
DCL::uetone(ux)

DCL::udsfmt("(f4.1)")
#ci = umax/ntone
ci = 0.2
DCL::udgcla(-ci*20,ci*20,ci)
DCL::uddclv(0)
=begin
[0.1,0.3,1,3].each{|x|
 DCL::udsclv(x,2,1,"#{x}",0.01)
 DCL::udsclv(-x,2,3,"-#{x}",0.01)
}
=end
DCL::udcntz(ux)

#DCL::usdaxs
DCL::uxaxdv("b",0.2,1)
DCL::uxaxdv("t",0.2,1)
DCL::uyaxdv("l",10,30)
DCL::uyaxdv("r",10,30)
#DCL::uxsttl("t","zonal-mean uhcos#{DCL::csgi(172)}",0)
DCL::uxsttl("b","time [s]",0)
DCL::uysttl("l","latitude",0)

vname = {"u-xm"=>"U","ucos-xm"=>"Ucos#{DCL::csgi(172)}","uhcos-xm"=>"Uhcos#{DCL::csgi(172)}"}[var]
DCL::sgtxzr(0.02,0.70,"zonal-mean #{vname}",0.02,0,-1,2)
if g
  DCL::sgtxzr(0.98,0.035,"Omega=#{omega.to_i}, gH=#{(g*h).to_i}, U=#{u0}, a=#{a.to_i}",0.012,0,1,1)
else
  DCL::sgtxzr(0.98,0.035,"Omega=#{omega.to_i}, gH=inf, U=#{u0}, a=#{a.to_i}",0.012,0,1,1)
end
#DCL::sgtxzr(0.98,0.015,"Contour interval=#{ci}",0.012,0,1,1)

DCL::grcls