# -*- coding: utf-8 -*-
require "numru/gphys"
include NumRu

#
# = 概要
#
# * 雨水量 (雨水混合比 x 密度) の鉛直積分: 
#
#     \int_z ( q_r * \rho ) dz
#
#   を計算する (ただし水平平均). 
#
# # たぶん時間平均する
#
# * チェック項目
#   * 用いるデータのファイル名
#   * 積分時間(tmax), タイムステップ数(tn)
# * 計算結果は thermal-moist_H2O-l-RainInt.nc に出力される
# * 最終的に出てくる量の変数を修正して出力するようにいるので, 
#   netCDF への出力がめんどくさい感じになっていると思う
#   # が, ひとまずこのままで
# 
# = めも
#
# * 新規作成: 2012/12/20
# * 最終更新: 2013/01/20
#

z1 = 0
z2 = ARGV[0].to_i
dz = 300
t2 = 4320000
#t2 = 2592000
dt = 18000
tn = t2 / tn

gphys1 = GPhys::IO.open("thermal-moist_H2O-l-RainAll.nc","H2O-l-RainAll")
gphys2 = GPhys::IO.open("thermal-moist_restart.nc","DensBZ")
gphys2 = gphys2.cut( 'x'=>0..256000 ).cut( "z"=>z1..z2 )

print "\n gphys1 * gphys2 :\n"
gphys = gphys1 * gphys2
p gphys, gphys.val

gphys = gphys.mean( 'y' )
gphys = gphys.mean( 'x' )
print "\n gpphys.mean(x).mean(y):\n"
p gphys
p ( gphys.data.get_att('units') )    # 次元チェック

print "\n gphys.sum(z) :\n"
gphys= gphys.sum( 'z' ) * dz
p gphys, gphys.val
p ( gphys.data.get_att('units') )    # 次元チェック

#-- NetCDF ファイルへの書き込み --#

t_f = gphys.val

  #-- 軸の設定 --#
t_a = VArray.new( NArray.sfloat(tn+1).indgen(0,t2/tn),
                    {"long_name"=>"time","units"=>"sec"},
                    "t" )
t = Axis.new.set_pos(t_a)

  #-- データ用の配列を準備 --#
data = VArray.new( NArray.sfloat(tn+1),
                   {"long_name"=>"Integral of H2O-l-Rain in z", "units"=>"Kg.m-2"},
                   "H2O-l-RainInt" )
t_f_gphys = GPhys.new( Grid.new(t), data )

  #-- 値を配列に入れる --#
t_f_gphys[0..tn] = t_f[0..tn]

  #-- netCDF ファイルの生成とファイルへの書き出し --#
outfile = NetCDF.create("thermal-moist_H2O-l-RainInt.nc")
GPhys::NetCDF_IO.write(outfile,t_f_gphys)

  #-- ファイルクローズ --#
outfile.close