#! /usr/bin/ruby
# -*- coding: utf-8 -*-
#
# 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 '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],
  ['--area',                      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'
outncfile = 'VerticalIntegrated' + $OPT_var + '.nc'

wsn = ($OPT_wsn||4)

area = ($OPT_area||'all')

  ## データ読み込み, 加工
  gv = gpopen dir + variable + ".nc"
  gps = gpopen dir + "Ps.nc"
  sigm = gpopen dir + "V.nc","sigm"

  return if gv.nil? || gps.nil?

  # 座標データの取得
  lon = gv.axis("lon")
  lat = gv.axis("lat")
  sig = gv.axis("sig")

  na_sig = sig.val
  na_lon = lon.val
  na_lat = lat.val
  na_sig_weight = GPhys::IO.open(ncfile, 'sig_weight').val
  na_lat_weight = GPhys::IO.open(ncfile, 'lat_weight').val

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

  data_name = 'IntegratedData'
  ofile = NetCDF.create(dir + data_name + '.nc')
  GPhys::NetCDF_IO.each_along_dims_write([gv,gps], ofile, 'time') { 
      |gdata,ps|  

      # 鉛直積算の計算
      gphys = gdata *  na_sig_weight.reshape(1, 1 ,nsig, 1) / na_sig_weight.sum
      vertical_integrated_data = gphys.sum('sig') * gps / Grav
      vertical_integrated_data = vertical_integrated_data

      # 空間平均の計算
      global_integrated_data = vertical_integrated_data * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum

      if (area == 'all') then
        data = global_integrated_data.mean('lon').sum('lat')
        data = data * 4.0 * Pi * RPlanet * RPlanet
      elsif (area == 'east') then
        data = global_integrated_data.cut('lon'=>180.0..360.0).mean('lon').sum('lat')
        data = data * 2.0 * Pi * RPlanet * RPlanet
      elsif (area == 'west') then
        data = global_integrated_data.cut('lon'=>0.0..179.9).mean('lon').sum('lat')
        data = data * 2.0 * Pi * RPlanet * RPlanet
      else
        p 'area is not set properly: area=', area
        exit
      end

      data = data.convert_units("kg")



      gvinteg = calc_msf(vwind,ps,sigm)


      [gvinteg]
      }
  ofile.close
  print "[#{data_name}](#{dir}) is created\n"




# ファイル出力
#outfile = NetCDF.create(outncfile)
#GPhys::NetCDF_IO.write( outfile, data )
#outfile.close

# (2013-11-16 石渡)
# 場合によっては
# [BUG] Segmentation fault
# になってしまって, やむなく gpview を使うこともあった.

exit



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( data, true, opt )

DCL.grcls