#
#  湿潤静的エネルギーをだす
#
##########################################################################

require "numru/gphys" # Gphys を使う
include NumRu

#-- 軸を読み込む/たぶんGphysオブジェクトを維持
#z = gphys.coord(2).val

#-- 定数

Grav = 9.8
CpDry = 1039.6423436145739742642
T_0   = 273.16 #[K], 水の融点

#--領域サイズ
#----大文字変数は定数(ruby)

#-- 幅設定

DelX = 500
DelZ =250
DelTime = 150
NT = (21600/150).to_i # どれだけ時間ステップか(NetCDFの軸を生成するために)
NX = 256
NZ = 120

#-- NArray.を用意

zTemp = NArray.sfloat(NX)

#x_fNl = NArray.sfloat(NX)
#x_fAl = NArray.sfloat(NX)

#-- NArray オブジェクト
#-- indgen(0,500)<- 0 から始めて 500 ずつ値を入れる
x_a = VArray.new(NArray.sfloat(NX).indgen(0,500), 
         {"long_name"=>"x-coodinate", "units"=>"m"}, # ハッシュ(ruby)
         "x")


#---- GPhys の軸のオブジェクトに格上げ
x = Axis.new.set_pos(x_a)

z_a = VArray.new(NArray.sfloat(NZ).indgen(0,250), 
         {"long_name"=>"z-coodinate", "units"=>"m"}, # ハッシュ(ruby)
         "z")

#---- GPhys の軸のオブジェクトに格上げ
z = Axis.new.set_pos(z_a)

#-- 時間のオブジェクト
t_a = VArray.new(NArray.sfloat(NT+1).indgen(0,150),
         {"long_name"=>"Time", "units"=>"sec"}, # ハッシュ(ruby)
          "t" )

#---- 軸のオブジェクトに格上げ
t = Axis.new.set_pos(t_a)

#-- 物理場オブジェクト
h_a = VArray.new(NArray.sfloat(NX,NZ,NT+1),
         {"long_name"=>"moist static energy", "units"=>"(1)"},
         "h")

#---- 軸のオブジェクトに格上げ
h_gphys = GPhys.new(Grid.new(x,z,t), h_a)

#-- netCDF づくり
ncFile = NetCDF.create("MSE.nc")
GPhys::NetCDF_IO.write_grid(ncFile, Grid.new(x,z,t))

###--- 必要な netCDF ファイルを open する ---###

gphysPTemp = GPhys::IO.open('./thermal_PTempAll.nc','PTempAll')
gphysExner = GPhys::IO.open('./thermal_ExnerAll.nc','ExnerAll')
gphysH2Og = GPhys::IO.open('./thermal_H2O-gAll.nc','H2O-gAll')

#-- 軸を読み込む/たぶんGphysオブジェクトを維持
tarray = gphysPTemp.coord(3).val
zarray = gphysPTemp.coord(2).val

### y
gphysPTemp = gphysPTemp.mean( 'y' )
gphysExner = gphysExner.mean( 'y' )
gphysH2Og = gphysH2Og.mean( 'y' )

for t in tarray
#t=0
p t 

### 時間でカット
   tPTemp = gphysPTemp.cut( 't'=>t )
   tExner = gphysExner.cut( 't'=>t )
   tH2Og = gphysH2Og.cut( 't'=>t )

   for z in zarray

### zでカット
      zPTemp = tPTemp.cut( 'z'=>z )
      zExner = tExner.cut( 'z'=>z )
      zH2Og = tH2Og.cut( 'z'=>z )
      zPTemp.units = '1'

### exp を計算するために NArray オブジェクトにする
      zPTemp = zPTemp.val
      zExner = zExner.val
      zH2Og = zH2Og.val

#-- 温度を出す

      zTemp = zPTemp*zExner

#-- 湿潤静止エネルギーh(T,q,z) = c_p T + L q + g z (L: 潜熱, q: 水蒸気混合比)
#-- 潜熱の計算
# [J/Kg]水の蒸発の潜熱. CReSSマニュアルからとってみた, Tは温度
# 簡易版 2.5x10^6, deepconv の値をつかってみてもいい
      latent  = 2.50078e6*(T_0/zTemp)**(0.167+(3.67*10e-4)*zTemp)
#      latent =  2.5e6

      h = (CpDry*zTemp) \
          + (latent * zH2Og) \
          + (Grav * z)
      k = ((z-125)/DelZ).to_i
      n = (t/DelTime).to_i
      h_gphys[0..NX-1,k,n] = h[0..NX-1]
   end
end

  h_gphys.data.rename!("h")
  h_gphys.data.set_att('long_name', "moist static energy")
  h_gphys.units = 'J/Kg'
  GPhys::NetCDF_IO.write(ncFile, h_gphys.rename('h'))  # netCDFファイルにライト
#rename で名前を変える
ncFile.close