#! /usr/bin/ruby
# -*- coding: utf-8 -*-
#
# 2015-3-17 石渡      鉛直積算は make_vertinteg_v3.rb で計算
#                     するようにする.
# 2015/02/19 石渡正樹 allocate memory error が出ないように改変.
#                     ファイルを出力してしまっているが, 
#                     ファイルを出力せずに図を描くようにしたい.
# 2014/10/03 石渡正樹
#
# = fig_globalinteg.rb
# * 全球積分値を求める.
#
# == USAGE
#   % fig_globalinteg.rb --var=H2OLiq --time_start=731 --time_end=1095 --wsn=2
# == Calculation Method
#  \int q d \sigma = \int q d dp/ps
#                  = 1/ps \int q \rho g (dz)  
#                  = g/ps \int \rho q (dz) 
# よって, 求めたいのは
# \int \rho q (dz)  = ps/g \int q d \sigma


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

ITR=1

## 定数設定
Grav = 9.8
RPlanet = 6.371e6
Pi = 3.1415926535897932

## オプション解析
parser = GetoptLong.new
parser.set_options(
       ###    global option   ###
       ['--var',                      GetoptLong::REQUIRED_ARGUMENT],
       ['--time_start',                      GetoptLong::REQUIRED_ARGUMENT],
       ['--time_end',                      GetoptLong::REQUIRED_ARGUMENT],
       ['--area',                      GetoptLong::REQUIRED_ARGUMENT],
       ['--wsn',                      GetoptLong::REQUIRED_ARGUMENT]
       )
parser.each_option do |name, arg|
  eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'"  # strage option value to $OPT_val
end

if $OPT_var == nil then
  p '--var must be specified.'
  exit
end
variable = $OPT_var
ncfile = $OPT_var + '.nc'

wsn = ($OPT_wsn||1)

time_start = ($OPT_time_start||1).to_f
time_end = ($OPT_time_end||9999).to_f

area = ($OPT_area||'all')

# データ入力
dir = "./"

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

sig = gpopen dir + ncfile,"sig"
lat = gpopen dir + ncfile,"lat"
lon = gpopen dir + ncfile,"lon"

sig_weight = gpopen dir + ncfile,"sig_weight"
lat_weight = gpopen dir + ncfile,"lat_weight"

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

na_lon = lon.val
na_lat = lat.val

na_lat_weight = lat_weight.val
na_sig_weight = sig_weight.val

nlat = na_lat.size

gp_fig = gqv.cut('time'=>time_start..time_end)

p gp_fig

# 緯度平均の計算
gp_fig = gp_fig * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum
gp_fig = gp_fig.sum('lat')
#gp_fig = gp_fig.cut('lon'=>lon_start..lon_end).mean('lon').sum('lat')

# 経度平均の計算
if (area == 'all') then
  gp_fig = gp_fig.mean('lon')
  gp_fig = gp_fig * 4.0 * Pi * RPlanet * RPlanet
elsif (area == 'east') then
  gp_fig = gp_fig.cut('lon'=>180.0..360.0).mean('lon')
  gp_fig = gp_fig * 2.0 * Pi * RPlanet * RPlanet
elsif (area == 'west') then
  gp_fig = gp_fig.cut('lon'=>0.0..180.0).mean('lon')
  gp_fig = gp_fig * 2.0 * Pi * RPlanet * RPlanet
else
  p 'area is not set properly: area=', area
  exit
end

# 鉛直平均
gp_fig = gp_fig * na_sig_weight


# 描画
DCL.gropn(wsn)

vx0,vx1,vy0,vy1 = 0.15,0.85,0.2,0.55
GGraph.set_fig 'itr'=>ITR, 'viewport'=>[vx0,vx1,vy0,vy1]

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

#  DCL.sgpset('lfull', true)     # 全画面表示
#  DCL.sgpset('lcntl', true)

  # オプションのデフォルト値の設定
  opt = {
    'legend'=>false, 'annotate'=>false,
    'title'=>'cloud water'
  }

   GGraph.line( gp_fig, true, opt )

DCL.grcls