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

#
# = 概要
#
# * 圧力方程式における熱膨張項を計算する
#   計算結果は thermal-moist_.nc に出力される
#   * 分子量効果の term のみ
#

#gphys = GPhys::NetCDF_IO.open('T.jan.nc', 'T')
#gphysPTempBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'PTempBZ')
gphysDensBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'DensBZ')
gphysPTempBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'PTempBZ')
gphysDensBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'DensBZ')
gphysVelSoudBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'VelSoundBZ')
gphysH2OgCond = GPhys::NetCDF_IO.open('thermal-moist_H2O-g_Cond.nc', 'H2O-g_Cond')
gphysH2OgSfcMassFlux = GPhys::NetCDF_IO.open('thermal-moist_H2O-g_SfcMassFlux.nc', 'H2O-g_SfcMassFlux')
gphysH2OlRainFall = GPhys::NetCDF_IO.open('thermal-moist_H2O-l-Rain_Fall.nc', 'H2O-l-Rain_Fall')
gphysH2OgAll = GPhys::NetCDF_IO.open('thermal-moist_H2O-gAll.nc', 'H2O-gAll')

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

xmin = 0
xmax = 256000
dx   = 1000
xn   = xmax / dx
zmin = 0
#zmax = 24000
#zmax = 30000
#zmax = 36000
zmax = 48000
dz   = 300
zn   = zmax / dz
tmin = 0
tmax = 86400 * 50
dt   = 1800
tn   = tmax / dt

CpDry = 1039.642343614574

#-- cut --#

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

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

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

gphysH2OgCond = gphysH2OgCond.mean( 'y' )
print "\n Condensation term of Water Vapor Mixing Ratio :\n"
p gphysH2OgCond.val

gphysH2OgSfcMassFlux = gphysH2OgSfcMassFlux.mean( 'y' )
print "\n Surface Mass Flux of Water Vapor Mixing Ratio :\n"
p gphysH2OgSfcMassFlux.val

gphysH2OlRainFall = gphysH2OlRainFall.mean( 'y' )
print "\n Fall Rain Term :\n"
p gphysH2OlRainFall.val

gphysH2OgAll = gphysH2OgAll.mean( 'y' )
print "\n Water Vaport Mixing Ratio :\n"
p gphysH2OgAll.val


switch = "n"

#if switch == "y"

  #-- 地表面フラックスを三次元の gphys オブジェクトにする --#

  xzt_f = gphysH2OgSfcMassFlux.val
  p xzt_f

    #-- 定数の設定 --#
  xmin = 500
  zmin = 150

    #-- 軸の設定 --#
  x_a = gphysH2OgCond.coord(0)
  x = Axis.new.set_pos(x_a)
  z_a = gphysH2OgCond.coord(1)
  z = Axis.new.set_pos(z_a)
  t_a = gphysH2OgCond.coord(2)
  t = Axis.new.set_pos(t_a)

  data = VArray.new( NArray.sfloat(xn,zn,tn+1),
                    {"long_name"=>"H2O-g surface mass flux", "units"=>"kg.m-2.s-1"},
                     "H2O-g_SfcMassFlux" )
  xzt_f_gphys = GPhys.new( Grid.new(x,z,t), data )
#  p "----------------------------"
  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,0..tn] = xzt_f[0..xn-1,0..tn]

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

#  p xzt_f_gphys

  # outfile を閉じる
  outfile.close
#----------------------------

#end

gphysH2OgSfcMassFlux = GPhys::NetCDF_IO.open('thermal-moist_2D-H2O-g_SfcMassFlux.nc', 'H2O-g_SfcMassFlux')

gphysH2OgSfcMassFlux = gphysH2OgSfcMassFlux / ( gphysDensBZ * dz )
print "\n Surface Water Vaport Mixing Ratio Flux :\n"
p gphysH2OgSfcMassFlux.val

#p "############################################"
#
##-- H2O-gAll の各高度での時間変化を計算 --#
#
####### この辺でエラーでる #####
#
#  #-- 配列を用意しておく --#
#
#aryH2OgAll = NArray.sfloat(xn,zn,tn+1)
#arydH2OgAlldt = NArray.sfloat(xn,zn,tn+1)
#
#aryH2OgAll = gphysH2OgAll.val
#
##for timestep in 0..tn
#for timestep in 0..2
#
#  t1 = timestep
#  t2 = timestep + 1 
#  aryH2OgAlldt = aryH2OgAll[0..xn-1,0..zn-1,t2] - aryH2OgAll[0..xn-1,0..zn-1,t1]  # 二次元配列に入れられる... 
#  p aryH2OgAlldt
#
##  p gphysH2OgAll.cut( 't'=>t2 ).val
##  xzt_f[0..xn-1,0..zn-1,timestep] = gphysH2OgAll.cut( 't'=>t2 ).val - gphysH2OgAll.cut( 't'=>t1 ).val - gphysH2OgCond.cut( 't'=>t1 ).val
##  gphysH2OgAll = gphysH2OgAll.cut( 't'=>t2 ) - gphysH2OgAll.cut( 't'=>t1 )
##  p gphysH2OgAll.val
##  xzt_f[0..xn-1,0..zn-1] = aryH2OgAll[0..xn-1,0..zn-1]
##  for i in 0..xn-1
##    for j in 0..zn-1
##       p aryH2OgAll[i,j]
###      xzt_f[i,j,timestep] = aryH2OgAll[i,j]
##    end
##  end
#
#end
##xzt_f[0..xn,0..zn,tn] = xzt_f[0..xn,0..zn,tn-1]
##
##p xzt_f
#
#exit
#
#data = VArray.new( NArray.sfloat(xn,zn,tn+1),
#                  {"long_name"=>"Tendency of H2O-gAll", "units"=>"kg.kg-1.s-1"},
#                   "dH2O-gAlldt" )
#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]
#xzt_f_gphys[0..xn-1,0..zn,0..tn] = xzt_f[0..xn-1,0..tn]
#
#    #-- netCDF ファイルの生成とファイルへの書き出し --#
##outfile = NetCDF.create("thermal-moist_ExnerThermExp.nc")
#outfile = NetCDF.create("thermal-moist_dH2O-gAlldt.nc")
#GPhys::NetCDF_IO.write(outfile,xzt_f_gphys)
#
#p xzt_f_gphys
#
## outfile を閉じる
#outfile.close
#

#-- 非断熱加熱の項を計算 (圧力, 移流・拡散なし) --#

gphysdpidt = ( gphysVelSoudBZ**2 / (CpDry * gphysDensBZ * gphysPTempBZ) )  * ( 0.609 * ( gphysH2OgCond + gphysH2OgSfcMassFlux ) - gphysH2OgCond - gphysH2OlRainFall )

print "\n Sumention of dqdt term in Pressure Eq. :\n"
p gphysdpidt.val
p gphysdpidt

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

  #-- 定数の設定 --#
xmin = 500
zmin = 150

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

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

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

t_a = gphysH2OgCond.coord(2)
#t_a = gphysPTempCond.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"},
                   "ExnerThermExpTerm" )
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("thermal-moist_Exner-ThermExpTerm-DqDt.nc")
GPhys::NetCDF_IO.write(outfile,xzt_f_gphys)

# outfile を閉じる
outfile.close