#
#  液体水静的エネルギーをだす
#   s_L = C_p T + gz - L q_L
##########################################################################

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

#-- 定数

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

#-- netCDF づくり
ncFile = NetCDF.create("LWSE.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')
gphysH2OlC = GPhys::IO.open('./thermal1_H2O-l-Cloud.nc','H2O-l-Cloud')
gphysH2OlR = GPhys::IO.open('./thermal1_H2O-l-Rain.nc','H2O-l-Rain')

z = gphysPTemp.axis(2).to_gphys

   yzPTemp = gphysPTemp
   yzExner = gphysExner
   yzH2Og = gphysH2Og
   yzH2OlC = gphysH2OlC
   yzH2OlR = gphysH2OlR

#-- 温度を出す

         yzTemp = yzPTemp*yzExner
         qL = yzH2OlC + yzH2OlR
#         qL = yzH2OlC    #-- 雲水だけバージョン

#-- (L: 潜熱, q: 水蒸気混合比)
#-- s_L = C_p T + gz - L q_L
#-- 潜熱の計算
# [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 * qL) \
             + (Grav * z)
         h_gphys = h

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