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

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

#-- 定数

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

#-- netCDF づくり
ncFile = NetCDF.create("MSE.nc")

###--- 必要な 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')

z = gphysPTemp.axis(2).to_gphys

   yzPTemp = gphysPTemp
   yzExner = gphysExner
   yzH2Og = gphysH2Og

#-- 温度を出す

         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)
         h_gphys = h

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