#!/usr/bin/env ruby
=begin

= gtstrm.rb : 南北風のデータから質量流線関数のデータを作成する.
  
= USAGE :

  $ gtstrm.rb -[hHo] <--output=[outputname]> [ncfile1] [ncfile2]

= note

* 同じ格子構造をもっているファイルでなくてはならない.

* 与えるデータは時間軸平均をかけたものを利用. (eddy は含まれていない.)

* 変数は一致しなくても演算は可能. 

* 結果として得られる変数は "strm".

* 流線関数値導出の式は 以下の式を導入. (TeX の書式で書いてます.)

  高度$z$ と南北方向$y$ および帯状平均した風速$[*]$を用いて, 連続の式を以下のように定義.
  \begin{align}
    \DP{[v]}{y} + \DP{[w]}{z} = 0
  \end{align}
  これに, 静水圧平衡の式
  \begin{align}
    \Dd p = \rho g \Dd z
  \end{align}
  を, 代入して z を圧力 p で書き表す. これを球面に拡張して, 以下を得る.
  \begin{align}
    \frac{\partial [\bar{v}] \cos \phi }{R\cos\phi\partial \phi} + 
             \DP{[\bar{\omega}]}{p} = 0. 
  \end{align}
  これを積分して
  \begin{align}
    \psi - \psi|_{p=0} & = \psi \\
                       & = \int_0^{1000}\Dd p \frac{2\pi R \cos\phi}{g}[\bar{v}]
  \end{align}
  本スクリプトでは, 上記方程式中の積分を台形公式にて差分化し, 解いた.

=end
## 保存する変数名
$new_vname = "strm"

require "getopts"
require "numru/ggraph"
include NumRu

## 追加メソッド定義ファイル
require "libgphys-n.rb"
## 属性設定ファイル
require "NCEP-ncattr.conf.rb"

#############  以下メインルーチン  ##############


# オプション解析 ----------------------------------------------------
unless getopts("hH", "help", "each", "o:", "output:", "c")
  print "#{$0}:illegal options. please exec this program with -h/--help. \n"
  exit 1
end


if ($OPT_h) || ($OPT_H) || ($OPT_help) || ARGV.size == 0 then
  print <<EOF

====================================================
                #{File.basename($0)} ver.0.1
====================================================

Describe:

  大気の南北風(v)のデータから, 子午面流線関数の値を求める.

  * データは 1 次元目が 緯度, 2 次元目が圧力でなければならない.

  * データフォーマットは netCDF に限る. 

  * 出力は, 与えられたデータと同様の格子構造をもつ netCDF フ
    ァイルである.

USAGE: #{File.basename($0)} <opt> [ncfile [ncfile ..] ] [var]
  
  * 引数解説

    [ncfile] : 対象とする netCDF ファイル名を指定. 
               (ex. mean-V_2001_01.nc)
               複数指定可能. その場合, 変数値を各格子毎に平均
               したものを出力する.

    [var]    : 南北風の変数名. できれば v であることが望ましい. 

  * オプション

    -h, -H, --help : このメッセージを表示.
    -o, --output   : 出力する netCDF ファイル名を指定. 
                     デフォルトは, "gtstrm.nc"
 
PARAMETER: 

  以下, コード中で与えているパラメータ.

  * r  :: 地球半径 => 6.370e6 (m)
  * g  :: 重力定数 => 9.8  (m/s^2)
  * pi :: 円周率   => #{Math::PI}

HISTORY:

  2003/12/06  created  by daktu32@ep.sci.hokudai.ac.jp

====================================================

EOF
exit 1
end



# 引数に与えた nc ファイルを格納する 配列.
file = Array.new

# 引数中の nc ファイルのみ格納
ARGV.each do |i| 
  if i =~ /.nc$/ then
    p i
    file << i
  end
end


# 変数 : デフォルトはヘッダ の一番最後の変数を自動的に選択
if ARGV[ARGV.size-1] =~ /.nc$/ then
  var = NetCDF.open(file[0], "r").var_names[-1]
else # 陽に指定も可能. その場合, 一番最後に与える.
  var = ARGV[ARGV.size-1]
end


# 南北風のデータじゃなかったらエラーを励起して終了.
# raise ArgumentError,"var-data is not nannbokuhuu." if var != "v"

# 保存ファイル名が与えられたら, その名前で保存. 
if ($OPT_o) then
  output = ($OPT_o).to_s
elsif ($OPT_output) then
  output = ($OPT_output).to_s
else 
  output = "gtstrm.nc"
end



# オプション解析 終了----------------------------------------------------

ask_overwrite(output)# 保存するファイルの上書きに関して問う. y or n で答える.

# メインルーチン

if ($OPT_each) # 複数ファイルを同時に扱うとき.
  file.each do |f|
    gphys = GPhys::NetCDF_IO.open(f, var)
    output = File.basename(f).sub("V_", "STRM_") if ($OPT_each)
    strm = gphys.mean(0).strm_function_ncep($new_vname)
    strm.save(output, $global_attr, var, $var_attr)          # nc ファイルに保存
  end
  exit
end

if file.size == 1 then  # 単体の nc ファイルの場合
  gphys = GPhys::NetCDF_IO.open(file[0].to_s, var)
  strm = gphys.mean(0).strm_function_ncep($new_vname) # 質量流線関数オブジェクト生成
  strm.save(output, $global_attr, $new_vname, $var_attr, $history)          # nc ファイルに保存
  exit
elsif ($OPT_c) 
  gphys = GPhys::NetCDF_IO.open(file, var)
else         
  vgarray = Array.new
  file.each do |f|
    vgarray << GPhys::NetCDF_IO.open(f.to_s, var)
  end
  vgphys = vgarray[0]
  1.upto(vgarray.length-1) {|n|
    vgphys += vgarray[n]
  }
  gphys = vgphys/file.length
end


p strm = gphys.mean(0).strm_function_ncep($new_vname)               # 質量流線関数オブジェクト生成
 strm.data.set_att('long_name', $var_attr['long_name'])
 strm.data.set_att('units', $var_attr['units'])

    if File.exist?(output)
      nc = NetCDF.open(output, 'a+')
    else
      nc = NetCDF.open(output, 'w')
    end
    GPhys::IO.write(nc, self)
    nc.close

    nc = NetCDF.open(output, "a+")
    nc.redef                      # define mode に移行

    # 大域属性付加
    $global_attr.to_a.each {|attr| 
      nc.put_att(attr[0], attr[1]) if attr[1]
    } if global_attr
    # history 付加
    now_hist = $global_attr["history"]
    if now_hist then
      hist = now_hist + "\n" + history
      nc.put_att("history", hist)
    else
      nc.put_att("history", history)
    end

    nc.close

    p "saved as #{output}"