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

# (2015-2-24 石渡) 長い時間でも計算できるようにした.
#                  でも擾乱成分 (時間平均からのずれ) の計算
#                  の実装はやってない.
#                  時間平均は時間ループの前に計算しておくと良い?

#DATA_HOME = "/home/noda/data"
#DATA_HOME = "/home/noda/data/exps_101220_diff_T21L16_sr_Omega_sweep"

require "/GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/dcpam.rb"
require "numru/ggraph"
include NumRu
require "pp"

wsn = 4

# ベクトルの描画間隔 (単位は格子点数)
xintv = 6
yintv = 2

#ITR=30   # 球
ITR=1

LON_CNT=125  # 球のときに使用

#xfact1 = 0.01
#yfact1 = 0.01
xfact1 = 1.0
yfact1 = 1.0

# 凡例の単位ベクトルの長さ
uxunit = 0.03
uyunit = 0.03

time_start = 1000
time_end   = 2000

#omegas = [0.0, 0.15, 0.5, 1.0]
#omegas = [1.0]
#omegas = [0.0]

omega = 1.0

type = 1  # 1: mean, 2: disturbance

cp = 1004.6

# 簡易
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 '--xfact1' then
    ARGV.shift
    xfact1 = ARGV[0]
    ARGV.shift
  when '--yfact1' then
    ARGV.shift
    yfact1 = ARGV[0]
    ARGV.shift
  when '--uxunit' then
    ARGV.shift
    uxunit = ARGV[0]
    ARGV.shift
  when '--uyunit' then
    ARGV.shift
    uyunit = ARGV[0]
    ARGV.shift
  when '--type' then
    ARGV.shift
    type = ARGV[0].to_i
    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
  else
    raise "ERROR: invalid option: #{ARGV[0]}"
  end
end



puts omega

if type == 1 then
  puts "all"
#elsif type == 2 then
#  puts "disturbance"
else
  raise "Type #{type} is not supported."
end

dir='./'

# geopotential 
temperature_file = dir + '/Temp.nc'
if File.exist?('Phi.nc') then
  puts("MESSAGE: Phi.nc exists. process skipped.")
else
  system("ruby /GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/geopotential_time.rb --var Temp.nc@Temp --range_time #{time_start}:#{time_end}")
end

# input
path = File.join(dir, 'U.nc')
na_sig_weight = GPhys::IO.open(path, 'sig_weight').val
nsig = na_sig_weight.size

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)

var = 'Temp'
path = File.join(dir, var+'.nc')
gp_h = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)*cp
gp_h = gp_h.convert_units('J kg-1')

# Phi
gp_phi = GPhys::IO.open('Phi.nc', 'Phi').cut('time'=>time_start..time_end)

gp_stengy = gp_h + gp_phi

sig = gpopen dir + "U.nc","sig"
na_sig = sig.val
nsig = na_sig.size

sig_weight = gpopen dir + "U.nc","sig_weight"
na_sig_weight = sig_weight.val

ofilex = NetCDF.create(dir + 'TrnsprtStEngyX' + '.nc')
ofiley = NetCDF.create(dir + 'TrnsprtStEngyY' + '.nc')

# 時間ループ
GPhys::NetCDF_IO.each_along_dims_write([gp_stengy,gp_u,gp_v], [ofilex,ofiley], 'time') {
  |gp_stengy,gp_u,gp_v|  

  gp_u_stengy = gp_u * gp_stengy
  gp_v_stengy = gp_v * gp_stengy

# # 擾乱成分
# if type == 2 then
#  gp_ud_qd = gp_u_q - gp_u * gp_q
#  gp_vd_qd = gp_v_q - gp_v * gp_q
# end

  # 鉛直積分
  gp_u_stengy = (gp_u_stengy * na_sig_weight.reshape(1,1,nsig)).sum('sig')
  gp_v_stengy = (gp_v_stengy * na_sig_weight.reshape(1,1,nsig)).sum('sig')

# 擾乱成分の場合の鉛直積分
# gp_ud_qd = (gp_ud_qd * na_sig_weight.reshape(1,1,nsig)).sum('sig')
# gp_vd_qd = (gp_vd_qd * na_sig_weight.reshape(1,1,nsig)).sum('sig')

  [gp_u_stengy,gp_v_stengy]
}
ofilex.close
ofiley.close

# 時間ループの後で時間平均が必要

sflux_x = GPhys::IO.open('TrnsprtStEngyX.nc', 'U').cut('time'=>time_start..time_end)
sflux_y = GPhys::IO.open('TrnsprtStEngyY.nc', 'V').cut('time'=>time_start..time_end)

sflux_x_tmean = sflux_x.mean('time')
sflux_y_tmean = sflux_y.mean('time')

###################################################

# 描画
DCL.gropn(wsn)

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

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

# 球の時
#GGraph.set_fig 'itr'=>ITR, 'viewport'=>[0.15,0.95,0.2,0.8], 'map_axis'=>[LON_CNT, 0, 0]
#GGraph.set_fig 'viewport'=>[0.15,0.95,0.2,0.8]

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


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

DCL.ugpset('XFACT1', xfact1)
DCL.ugpset('YFACT1', yfact1)

DCL.ugpset('UXUNIT', uxunit) # size of unit vector by dimentional value
DCL.ugpset('UYUNIT', uyunit)

DCL::ugrset('VXUNIT', 0.1)  # size of unit vector
DCL::ugrset('VYUNIT', 0.1)

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

DCL.ugpset('LMSG', false)   # no message

DCL.ugpset('INDEX', 3)

opt_vector = {
  'xintv'=>xintv, 'yintv'=>yintv, 'title'=>'',
  'flow_vect'=>false, 'annotate'=>false, 'unit_vect'=>true
}

# mean
GGraph.vector( sflux_x_tmean, sflux_y_tmean, true, opt_vector)

# 擾乱の場合
# GGraph.vector( gp_ud_qd, gp_vd_qd, true, opt_vector)

DCL.grcls