# -*- 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 = 256000
z1 = 0
z2 = ARGV[0].to_i
dz = 300
t2 = 4320000
#t2 = 2592000
dt = 1800
tn = t2 / dt

gphysH2OgAll = GPhys::IO.open("thermal-moist_H2O-gAll.nc","H2O-gAll")
gphysDensBZ  = GPhys::IO.open("thermal-moist_restart.nc","DensBZ")
gphysDensBZ  = gphysDensBZ.cut( 'x'=>x1..x2 ).cut( "z"=>z1..z2 )

print "\n gphys1 * gphys2 :\n"
gphysWetH2OgAll = gphysH2OgAll * gphysDensBZ
p gphysWetH2OgAll, gphysWetH2OgAll.val

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

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


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

t_f = gphysTotalWetH2OgAll.val

  #-- 軸の設定 --#
t_a = VArray.new( NArray.sfloat(tn+1).indgen(0,dt),
                    {"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("thermal-moist_H2O-gInt.nc")
GPhys::NetCDF_IO.write(outfile,t_f_gphys)

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