#! /usr/bin/ruby
# -*- coding: utf-8 -*-
#
# 2015-07-08 石渡     ソースの中身を整理
# 2015-03-16 石渡     出力ファイルの変数名を VertInt + var にした.
#                     軸データを入れるべし.
#
# 2015/02/19 石渡正樹 allocate memory error が出ないように改変.
#                     ファイルを出力してしまっているが, 
#                     ファイルを出力せずに図を描くようにしたい.
# 2014/10/03 石渡正樹
#
# = make_vertiteg.rb
# == 概要
# * 鉛直積算値を計算してその結果をファイルに出力する.
#   出力されるファイル名の形式は VertInt(変数名).nc 
#
# == 使用例
#   % fig_vertinteg.rb --var=H2OLiq --time_start=731 --time_end=1095
#
# == 鉛直積算値計算方法
# * 鉛直積算
#   \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 'getoptlong'
require "numru/ggraph"
include NumRu

## 定数設定
Grav = 9.8

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

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

dir = "./"

var = GPhys::IO.open dir+ncfile,$OPT_var
ps = GPhys::IO.open dir+'Ps.nc', 'Ps'
sigm = GPhys::IO.open dir+ncfile, 'sigm'

sig = GPhys::IO.open dir+ncfile, 'sig'
lat = GPhys::IO.open dir+ncfile, 'lat'
lon = GPhys::IO.open dir+ncfile, 'lon'

sig_weight = GPhys::IO.open dir+ncfile, 'sig_weight'
lat_weight = GPhys::IO.open dir+ncfile, 'lat_weight'
lon_weight = GPhys::IO.open dir+ncfile, 'lon_weight'

return if var.nil? || ps.nil?

na_sig_weight = sig_weight.val
nsig = sig.val.size

data_name = 'VertInt' + $OPT_var
ofile = NetCDF.create(dir + data_name + '.nc')

# 時間ループ
GPhys::NetCDF_IO.each_along_dims_write([var,ps], ofile, 'time') { 
  |var,ps|  
        
  time = var.axis("time")    
        
  vintvar = var.copy
  vintvar[false] = 0
        
  alph = var * ps / Grav 

  # 鉛直積算の計算
  vintvar = alph *  na_sig_weight.reshape(1, 1 ,nsig, 1)/na_sig_weight.sum
  data = vintvar.sum('sig')
  data.units = 'kg m-2'
  data.name = data_name
  [data]
}

# 軸データの書き出し
GPhys::NetCDF_IO.write( ofile, lon_weight )
GPhys::NetCDF_IO.write( ofile, lat_weight )
GPhys::NetCDF_IO.write( ofile, sig_weight )

ofile.close

exit