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

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

# TODO
# (2014-11-11 石渡) 潜熱エネルギーの計算部分をメソッドに変更する.

require "/GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/dcpam.rb"

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

Grav = 9.8

RPlanet = 6.371e6

Pi = 3.1415926535897932

Latent = 2500000
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

# 簡易
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 = File.join(DATA_HOME, 
#      "101220_diff-DryAir-Vapor_T21L16_sr_Omega" + omega.to_s + "_dt20m")

dir='./'


# 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 = 'QVap'
#path = File.join(dir, var+'.nc')
#gp_q = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)*Latent

# gpq
gqvap = gpopen dir + "QVap.nc"
gu = gpopen dir + "U.nc"
gv = gpopen dir + "V.nc"

ofilex = 'TrnsprtLqX.nc'
ofiley = 'TrnsprtLqY.nc'


#ここから

gqv = gpopen dir + "H2OLiq.nc"
gps = gpopen dir + "Ps.nc"
sigm = gpopen dir + "H2OLiq.nc","sigm"

sig = gpopen dir + "H2OLiq.nc","sig"
lat = gpopen dir + "H2OLiq.nc","lat"
lon = gpopen dir + "H2OLiq.nc","lon"

sig_weight = gpopen dir + "H2OLiq.nc","sig_weight"
lat_weight = gpopen dir + "H2OLiq.nc","lat_weight"

return if gqv.nil? || gps.nil?

na_sig = sig.val
na_lon = lon.val
na_lat = lat.val

na_sig_weight = sig_weight.val
na_lat_weight = lat_weight.val

nsig = na_sig.size
nlat = na_lat.size

ofilex = NetCDF.create(dir + 'TrnsprtLqX' + '.nc')
ofiley = NetCDF.create(dir + 'TrnsprtLqY' + '.nc')

# 時間ループ
GPhys::NetCDF_IO.each_along_dims_write([gu,gv,gqvap], [ofilex,ofiley], 'time') { 
  |gp_u,gp_v,qvap|  
        
  time = qvap.axis("time")    
        
  gp_u_q = gp_u * qvap
  gp_v_q = gp_v * qvap

#  # 擾乱成分
#  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_q = (gp_u_q * na_sig_weight.reshape(1,1,nsig)).sum('sig')
  gp_v_q = (gp_v_q * 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_q,gp_v_q]
}
ofilex.close
ofiley.close


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

# フラックス鉛直積分データの読み込み
#qflux_x = gpopen dir + "TrnsprtLqX.nc"
#qflux_y = gpopen dir + "TrnsprtLqY.nc"

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

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

qflux_x_tmean = qflux_x.mean('time')*Latent
qflux_y_tmean = qflux_y.mean('time')*Latent


#exit


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


# 描画
DCL.gropn(wsn)

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( qflux_x_tmean, qflux_y_tmean, true, opt_vector)

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

DCL.grcls