# -*- coding: euc-jp -*-
require "numru/gphys"
include NumRu

#
# = 概要
#
# * 圧力方程式における熱膨張項を計算する
#   計算結果は thermal-moist_Exner-ThermExpTerm-DthetaDt.nc に出力される
#   * 今のところは, 非断熱加熱の term だけ
#   * 地表面熱フラックス入れてない版
#

#gphys = GPhys::NetCDF_IO.open('T.jan.nc', 'T')
gphys1  = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'PTempBZ')
gphys2  = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'DensBZ')
gphys3  = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'VelSoundBZ')
#gphys11 = GPhys::NetCDF_IO.open('thermal-moist_PTempCond.nc', 'PTempCond')
#gphys12 = GPhys::NetCDF_IO.open('thermal-moist_PTempDisp.nc', 'PTempDisp')
#gphys13 = GPhys::NetCDF_IO.open('thermal-moist_PTempRad.nc', 'PTempRad')
#gphys14 = GPhys::NetCDF_IO.open('thermal-moist_PTempTurb.nc', 'PTempTurb')
gphys11 = GPhys::NetCDF_IO.open('thermal-moist_PTempCond_IntValue30day.nc', 'PTempCond')
gphys12 = GPhys::NetCDF_IO.open('thermal-moist_PTempDisp_IntValue30day.nc', 'PTempDisp')
gphys13 = GPhys::NetCDF_IO.open('thermal-moist_PTempRad_IntValue30day.nc', 'PTempRad')
gphys14 = GPhys::NetCDF_IO.open('thermal-moist_PTempTurb_IntValue30day.nc', 'PTempTurb')


#-- 必要な定数を設定 --#

xmin = 0
xmax = 256000
zmin = 0
#zmax = 24000
zmax = 30000
#zmax = 36000
#zmax = 48000
dz   = 300
tmax = 2592000

CpDry = 1039.642343614574

#-- cut --#

gphys1  = gphys1.cut( 'y'=>500 )
gphys1  = gphys1.cut( 'x'=>xmin..xmax )
gphys1  = gphys1.cut( 'z'=>zmin..zmax )
print "\n Basic field of Potential Temperature :\n"
p gphys1.val

gphys2  = gphys2.cut( 'y'=>500 )
gphys2  = gphys2.cut( 'x'=>xmin..xmax )
gphys2  = gphys2.cut( 'z'=>zmin..zmax )
print "\n Basic field of Density :\n"
p gphys2.val

gphys3  = gphys3.cut( 'y'=>500 )
gphys3  = gphys3.cut( 'x'=>xmin..xmax )
gphys3  = gphys3.cut( 'z'=>zmin..zmax )
print "\n Basic field of Sound Velocity :\n"
p gphys3.val

gphys11 = gphys11.cut( 'y'=>500 )
gphys11 = gphys11 * tmax
print "\n Sumention of Condensation term in Potential Temperature :\n"
p gphys11.val

gphys12 = gphys12.cut( 'y'=>500 )
gphys12 = gphys12 * tmax
print "\n Sumention of Dispation term in Potential Temperature :\n"
p gphys12.val

gphys13 = gphys13.cut( 'y'=>500 )
gphys13 = gphys13 * tmax
print "\n Sumention of Radiation term in Potential Temperature :\n"
p gphys13.val

gphys14 = gphys14.cut( 'y'=>500 )
gphys14 = gphys14 * tmax
print "\n Sumention of Turbulence term in Potential Temperature :\n"
p gphys14.val

#-- 非断熱加熱の項を計算 --#

gphysdtheta = ( gphys2 / gphys1 ) * ( gphys11 + gphys12 + gphys13 + gphys14 )

print "\n Sumention of diabatic heat term :\n"
p gphysdtheta.val

gphysdrho = -gphysdtheta 

gphysdpidt = ( gphys3**2 / (CpDry * gphys1) )  * ( 1 / gphys1 ) * ( gphys11 + gphys12 + gphys13 + gphys14 )
print "\n Sumention of diabatic heat term in Pressure Eq. :\n"
p gphysdpidt.val

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

  #-- 定数の設定 --#
xmin = 500
xn = 256
zmin = 150
#zn = 80
zn = 100
#zn = 120
#zn = 160
tmin = 0
#tn = 2400
#tn = 1440
tn = 1

  #-- 値を配列に入れておく --#
xzt_f = gphysdrho.val

  #-- 軸の設定 --#
x_a = gphys11.coord(0)
x = Axis.new.set_pos(x_a)

z_a = gphys11.coord(1)
z = Axis.new.set_pos(z_a)

t_a = gphys11.coord(2)
t = Axis.new.set_pos(t_a)

  #-- drho 用 --#
    #-- データ用の配列を準備 --#
data = VArray.new( NArray.sfloat(xn,zn,tn+1),
                  {"long_name"=>"Contribution of Thermal Expansion Term in Density", "units"=>"Kg.m-3"},
                   "DensThermExpTerm" )
xzt_f_gphys = GPhys.new( Grid.new(x,z,t), data )
p xzt_f_gphys

    #-- 値を配列に入れる --#
xzt_f_gphys[0..xn-1,0..zn-1,0..tn] = xzt_f[0..xn-1,0..zn-1,0..tn]
p xzt_f_gphys

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

# outfile を閉じる
outfile.close


  #-- dpidt 用 --#
    #-- 値を配列に入れておく --#
xzt_f = gphysdpidt.val
p xzt_f


    #-- データ用の配列を準備 --#
data = VArray.new( NArray.sfloat(xn,zn,tn+1),
                  {"long_name"=>"Contribution of Thermal Expansion Term in Pressure Eq.", "units"=>"1"},
                   "ExnerThermExpTerm" )
xzt_f_gphys = GPhys.new( Grid.new(x,z,t), data )
p xzt_f_gphys

    #-- 値を配列に入れる --#
xzt_f_gphys[0..xn-1,0..zn-1,0..tn] = xzt_f[0..xn-1,0..zn-1,0..tn]
p xzt_f_gphys

    #-- netCDF ファイルの生成とファイルへの書き出し --#
#outfile = NetCDF.create("thermal-moist_ExnerThermExp.nc")
outfile = NetCDF.create("thermal-moist_Exner-ThermExpTerm-DthetaDt.nc")
GPhys::NetCDF_IO.write(outfile,xzt_f_gphys)


# outfile を閉じる
outfile.close

exit

#-------------------------------------------------#

# 変数の rename など
#gphysHumidity.data.rename!("RevHumid")
#gphysHumidity.data.set_att('long_name', "rerative humidity")

# 出力らへん
#outfile = NetCDF.create('tmp.nc')
outfile = NetCDF.create('zzz_test-thermal-moist_NetsuboutyouTerm.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, ( gphysdrho ).rename('ThermExpTerm') )
#GPhys::NetCDF_IO.write( outfile, gphys )