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

#
# = 概要
#
# * 圧力方程式における熱膨張項を計算する
#   計算結果は zzz_test-thermal-moist_NetsuBoutyouTerm.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')
gphys10 = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll_IntValue30day.nc', 'ExnerAll')
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')
gphys15 = GPhys::NetCDF_IO.open('thermal-moist_SfcHeatFlux_IntValue30day.nc', 'SfcHeatFlux')


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

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 )
gphys1  = gphys1.cut( 'z'=>150 )
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 )
gphys2  = gphys2.cut( 'z'=>150 )
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 )
gphys3  = gphys3.cut( 'z'=>150 )
print "\n Basic field of Sound Velocity :\n"
p gphys3.val

gphys10 = gphys10.cut( 'y'=>500 )
gphys10 = gphys10.cut( 'z'=>150 )
print "\n Exner function :\n"
p gphys10.val

gphys11 = gphys11.cut( 'y'=>500 )
gphys11 = gphys11.cut( 'z'=>150 )
print "\n Condensation term of Potential Temperature :\n"
p gphys11.val

gphys12 = gphys12.cut( 'y'=>500 )
gphys12 = gphys12.cut( 'z'=>150 )
print "\n Dispation term of Potential Temperature :\n"
p gphys12.val

gphys13 = gphys13.cut( 'y'=>500 )
gphys13 = gphys13.cut( 'z'=>150 )
print "\n Radiation term of Potential Temperature :\n"
p gphys13.val

gphys14 = gphys14.cut( 'y'=>500 )
gphys14 = gphys14.cut( 'z'=>150 )
print "\n Turbulence term of Potential Temperature :\n"
p gphys14.val

gphys15 = gphys15.cut( 'y'=>500 )
gphys15 = gphys15 / (CpDry * gphys2 * gphys10 * dz )
print "\n Surface Potential Temperature Flux :\n"
p gphys15.val


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

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

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

gphysdrho = -gphysdtheta 

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

#-- 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_a = gphys11.coord(1)
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"},
#                   "DensThermExp" )
#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_DensThermExp.nc")
#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"},
#                   "ExnerThermExp" )
#xzt_f_gphys = GPhys.new( Grid.new(x,z,t), data )
data = VArray.new( NArray.sfloat(xn,tn+1),
                  {"long_name"=>"Contribution of Thermal Expansion Term in Pressure Eq.", "units"=>"1"},
                   "ExnerThermExp" )
xzt_f_gphys = GPhys.new( Grid.new(x,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]
xzt_f_gphys[0..xn-1,0..tn] = xzt_f[0..xn-1,0..tn]
p xzt_f_gphys

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

# outfile を閉じる
outfile.close