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

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

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

#-- 定数

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

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

#-- 幅設定

DelX = 400
DelY = 400
DelZ = 200
DelTime = 120
NT = (10800/120).to_i # どれだけ時間ステップか(NetCDFの軸を生成するために)
NX = 50
NZ = 50

#-- 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,DelX), 
         {"long_name"=>"x-coodinate", "units"=>"m"}, # ハッシュ(ruby)
         "x")

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

#-- y
y_a = VArray.new(NArray.sfloat(NX).indgen(0,DelY), 
         {"long_name"=>"y-coodinate", "units"=>"m"}, # ハッシュ(ruby)
         "y")

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

#-- z
z_a = VArray.new(NArray.sfloat(NZ).indgen(0,DelZ), 
         {"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,DelTime),
         {"long_name"=>"Time", "units"=>"sec"}, # ハッシュ(ruby)
          "t" )

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

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

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

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

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

gphysPTemp = GPhys::IO.open('./thermal1_PTempAll.nc','PTempAll')
gphysExner = GPhys::IO.open('./thermal1_ExnerAll.nc','ExnerAll')
gphysH2Og = GPhys::IO.open('./thermal1_H2O-gAll.nc','H2O-gAll')

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

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 )

### yでカット
      for y in yarray
         yzPTemp = zPTemp.cut( 'y'=>y )
         yzExner = zExner.cut( 'y'=>y )
         yzH2Og = zH2Og.cut( 'y'=>y )

### exp を計算するために NArray オブジェクトにする
         yzPTemp = yzPTemp.val
         yzExner = yzExner.val
         yzH2Og = yzH2Og.val

#-- 温度を出す

         yzTemp = yzPTemp*yzExner

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

         h = (CpDry*yzTemp) \
             + (latent * yzH2Og) \
             + (Grav * z)
         l = ((y-200)/DelY).to_i
         k = ((z-100)/DelZ).to_i
         n = (t/DelTime).to_i
         h_gphys[0..NX-1,l,k,n] = h[0..NX-1]
      end
   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)  # netCDFファイルにライト
#rename で名前を変える
ncFile.close