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

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

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

gphys1 = GPhys::IO.open("thermal1_H2O-g.nc","H2O-g")
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

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') )    # 次元チェック


#-- 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 Disturbance in z", "units"=>"Kg.m-2"},
                   "H2O-gDisInt" )
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-gDisInt.nc")
GPhys::NetCDF_IO.write(outfile,t_f_gphys)

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