#!/usr/bin/ruby
# -*- coding: utf-8 -*-

# 2014/09/14 Ishiwatari  modified

# Satoh et al (1995) の PV の式は
#  PV = - g * 1/(a^2 cos phi) \DP[\theta]{l}{phi} \DP{\theta}{p} 
# 温位面で微分をとらないといけない.

# GPhys の PV 計算モジュールを使ってみる.
# 温位面を指定する.

require "numru/ganalysis"
require "numru/ggraph"
include NumRu
require "pp"

wsn = 4

RPlanet = 6.371 * 10 ** 6  # [m]
Grav = 9.8   # [m/s2]

PI=3.141592

title = ""

# vapor の描画範囲
val_min = 0
val_max = 3e9
#val_interval = (val_max - val_min) / 30.0
val_interval = 1e8

#ITR=30   # 球
ITR=1   # 平面プロット

LON_CNT=125  # 球のときに使用

time_start = 400
time_end   = 500

omega = 7.292106590880652e-5

dir = nil

# 簡易
while ARGV.size > 0
  case ARGV[0]
  when '--omega' then
    ARGV.shift
    omega = ARGV[0]
    ARGV.shift
  when '--wsn' then
    ARGV.shift
    wsn = ARGV[0]
    ARGV.shift
  when '--xintv' then
    ARGV.shift
    xintv = ARGV[0]
    ARGV.shift
  when '--yintv' then
    ARGV.shift
    yintv = ARGV[0]
    ARGV.shift
  when '--range_time' then
    ARGV.shift
    tmp = ARGV[0].split(':')
    time_start = tmp[0].to_f
    time_end   = tmp[1].to_f
    ARGV.shift
  when '--range_val' then
    ARGV.shift
    tmp = ARGV[0].split(':')
    val_min = tmp[0].to_f
    val_max = tmp[1].to_f
    ARGV.shift
  when '--cint' then
    ARGV.shift
    val_interval = ARGV[0]
    ARGV.shift
  else
    raise "ERROR: invalid option: #{ARGV[0]}"
  end
end

dir='.'

# データ読み込み
path = File.join(dir, 'U.nc')

var = 'Temp'
path = File.join(dir, var+'.nc')
gp_temp = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)

var = 'U'
path = File.join(dir, var+'.nc')
gp_u = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)

var = 'V'
path = File.join(dir, var+'.nc')
gp_v = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)

lon  = GPhys::IO.open(path, 'lon')
lat  = GPhys::IO.open(path, 'lat')
sig  = GPhys::IO.open(path, 'sig')
time = GPhys::IO.open(path, "time").cut('time'=>time_start..time_end)

nlon  = lon.val.size
nlat  = lat.val.size
nsig  = sig.val.size
ntime = time.val.size

var = 'Ps'
path = File.join(dir, var+'.nc')
gp_ps = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)

na_gp_ps = gp_ps.val
na_sig = sig.val

pressure = NArray.sfloat(nlon,nlat,nsig,ntime).fill(0.0)

for i in 0..nlon-1
  for j in 0..nlat-1
    for k in 0..nsig-1
      for t in 0..ntime-1
        pressure[i,j,k,t] = na_gp_ps[i,j,t]*na_sig[k]
      end
    end
  end
end

gp_press = VArray.new( pressure,
                   {"long_name"=>"Pressure", "units"=>"Pa"},
                   "Pressure" )

theta = GAnalysis::Met.temp2theta(gp_temp,prs=gp_press)

p theta


# 角運動量の計算
lat_rad = lat*PI/180.0

cos_lat_rad = lat_rad.cos
na_cos_lat_rad = cos_lat_rad.val

na_cos_lat_rad_4d = na_cos_lat_rad.reshape(1,nlat,1,1)

#p omega*RPlanet*RPlanet*na_cos_lat_rad_4d*na_cos_lat_rad_4d

gp_angularmom = RPlanet*gp_u*lat_rad.cos+omega*RPlanet*RPlanet*na_cos_lat_rad_4d*na_cos_lat_rad_4d

gp_angularmom_mean = gp_angularmom.mean('time').mean('lon')

# p gp_angularmom_mean

DCL.gropn(wsn)

DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7)

GGraph.set_axes('xtickint'=>10, 'xlabelint'=>30)
GGraph.set_axes('ytickint'=>10, 'ylabelint'=>30)

# same as gpview
vx0,vx1,vy0,vy1 = 0.15,0.85,0.2,0.55
GGraph.set_fig 'itr'=>ITR, 'viewport'=>[vx0,vx1,vy0,vy1], 'map_axis'=>[LON_CNT, 0, 0]

DCL.uscset('cyspos', 'B' )              # move unit y axis

DCL.udlset('LMSG', false)

DCL::uxmttl('T', '', 0.0)  # title (large character)

#DCL.ugpset('LNRMAL', false)  # scale vector automatically
DCL.ugpset('LUMSG', false)   # no message

DCL.ugpset('INDEX', 9)

opt_contour = {
  'min'=>val_min, 'max'=>val_max, 'interval'=>val_interval,
  'title'=>title, 'annotate'=>false
}

GGraph.tone_and_contour( gp_angularmom_mean, true, opt_contour)
GGraph.color_bar('landscape'=>true)

DCL.grcls