# -*- coding: euc-jp -*-

# エクスナー関数と標準圧力から気圧を計算し, netCDF ファイルにはくスクリプト

require "numru/gphys"
include NumRu

gphysExner = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll.nc', 'ExnerAll')
gphysPressure0 = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'PressBZ')

x1 = 0
x2 = 256000
z1 = 0
z2 = 48000

CpDry    = 1039.64234361457397
GasRDry  = 296.785690317661306

gphysPressure0 = gphysPressure0.cut( 'x'=>x1..x2 )
gphysPressure0 = gphysPressure0.cut( 'z'=>z1..z2 )

# p = p_0 * \pi^{c_{pd} / R_d}
gphysPressure = gphysPressure0 * gphysExner ** (CpDry / GasRDry)

# 変数の rename など
gphysPressure.data.rename!("PressAll")
gphysPressure.data.set_att('long_name', "Pressure")

# 出力らへん
outfile = NetCDF.create('thermal-moist_PressAll.nc')
GPhys::NetCDF_IO.write( outfile, gphysPressure )

# outfile を閉じる
outfile.close

exit


##### thermal-moist_TempAll_rank000001.nc を作成 #####

#gphys = GPhys::NetCDF_IO.open('T.jan.nc', 'T')
gphys1 = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll_rank000001.nc', 'ExnerAll')
gphys2 = GPhys::NetCDF_IO.open('thermal-moist_PTempAll_rank000001.nc', 'PTempAll')

# 温度 = エクスナー関数 x 温位
gphys = gphys1 * gphys2

# 変数の rename など
gphys.data.rename!("TempAll")
gphys.data.set_att('long_name', "Temperature")

# 出力らへん
#outfile = NetCDF.create('tmp.nc')
outfile = NetCDF.create('thermal-moist_TempAll_rank000001.nc')
#GPhys::NetCDF_IO.write( outfile, gphys.cut('level'=>1000..250).mean(0) )
#GPhys::NetCDF_IO.write( outfile, gphys.cut('level'=>1000..250).mean(0,1).rename('T00') )
#GPhys::NetCDF_IO.write( outfile, ( gphys1*gphys2 ).rename('Temp') )
GPhys::NetCDF_IO.write( outfile, gphys )

# outfile を閉じる
outfile.close