#! /usr/bin/ruby
# -*- coding: utf-8 -*-
#
# 2007/09/28 石渡正樹
# 2013/05/26 石渡正樹
#
# = fig_vertinteg-cloudwater.rb
# * 鉛直積算雲水量の図を作る.
#
# == USAGE
#   % fig_vertinteg-cloudwater.rb --time_start=731 --time_end=1096 --wsn=2

require 'getoptlong'
require "numru/ggraph"
include NumRu

## 定数設定
Grav = 9.8

## オプション解析
parser = GetoptLong.new
parser.set_options(
                   ###    global option   ###
                   ['--time_start',                      GetoptLong::REQUIRED_ARGUMENT],
                   ['--time_end',                      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

wsn = ($OPT_wsn||4)

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

## データ読み込み, 加工
H2OLiq = GPhys::IO.open('H2OLiq.nc', 'H2OLiq').cut('time'=>time_start..time_end)
Ps = GPhys::IO.open('Ps.nc', 'Ps').cut('time'=>time_start..time_end)

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

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

na_sig_weight = GPhys::IO.open('H2OLiq.nc', 'sig_weight').val

nsig = na_sig.size
nlon = na_lon.size
nlat = na_lat.size
ntime = na_time.size

gphys = H2OLiq *  na_sig_weight.reshape(1, 1 ,nsig, 1) / na_sig_weight.sum
integrated_water = gphys.sum('sig') / Ps * Grav
integrated_water = integrated_water.mean('time')

p integrated_water

outfile = NetCDF.create('IntegH2OLiq.nc')
GPhys::NetCDF_IO.write( outfile, integrated_water )
outfile.close



# (2013-11-16 石渡)
# 以下で図を描こうとすると, 
# [BUG] Segmentation fault
# になってしまう.
# やむなく gpview を使うことにした.



exit

DCL.gropn(wsn)

GGraph.set_fig( 'itr'=> 1, 'viewport'=>[0.2,0.8,0.3,0.6] )

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

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

#   GGraph.tone_and_contour( integrated_water, true )
#   GGraph.contour( integrated_water, true )
   GGraph.tone( integrated_water, true )

DCL.grcls


