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

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

x1 = 0
x2 = 20000
y1 = 0
y2 = 20000
z1 = 0
z2 = 20000

#gp = GPhys::IO.open "air.2012-01.nc","air"
#gp = gp.cut("level"=>925..20, "lat"=>40)[true,{0..-1=>2},true]
#     #^サブセット指定 (この時点ではまだ読み込まれない)
#     # 後ろの [ ] があると, gphys の依存変数 z が一格子点飛びに表示された
gphys1 = GPhys::IO.open("thermal1_H2O-gAll.nc","H2O-gAll")
gphys2 = GPhys::IO.open("thermal_restart.nc","DensBZ")
gphys2  = gphys2.cut( 'x'=>x1..x2 ).cut( "z"=>z1..z2 )

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

#print "gp :\n"
#p gp
#print "pressure levels:\n"
#p gp.coord("level").val
#print "\ngp.mean(0) :\n"
#p gp.mean(0)
#print "\ngp.mean('lon','time') :\n"
#p gp.mean('lon','time')
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( 'x' )
#gphys= gphys.sum( 'y' )
gphys= gphys.sum( 'z' )
p gphys, gphys.val
p ( gphys.data.get_att('units') )    # 次元チェック

#print "\nmean,max,min,stddev values:\n"
#p gp.mean('lon','time').val, gp.max('lon','time').val, 
#  gp.min('lon','time').val, gp.stddev('lon','time').val


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

#tmax = 4320000
tmax = 10800
#tn = 2400
tn = 90
t_f = gphys.val

  #-- 軸の設定 --#
t_a = VArray.new( NArray.sfloat(tn+1).indgen(0,tmax/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-g in z", "units"=>"Kg.m-2"},
                   "H2O-gInt" )
t_f_gphys = GPhys.new( Grid.new(t), data )

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

  #-- テスト出力 --#
#p "-----------------------------"
#p t_f
#p t_f_gphys

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

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