#
#  相当温位をだす
#
##########################################################################

require "numru/gphys" # Gphys を使う
include NumRu
include NMath # NArray で計算できるように


#-- 定数
Grav = 9.8
CpDry = 1004.0 # 標準的 CpDry
T_0   = 273.16 #[K], 水の融点
P_0   = 1010.0 #deepconv 計算中での基準気圧 hPa
R_d  = 287.0 # J/kg 乾燥気体定数
R_v   = 461.0 # J/kg 水蒸気の気体定数

#-- netCDF づくり
ncFile = NetCDF.create("EPT.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でカット
      yzPTemp = gphysPTemp
      yzExner = gphysExner
      yzH2Og = gphysH2Og

#-- 温度を計算

         yzTemp = yzPTemp*yzExner

#-- 圧力を計算
         p = P_0*(yzExner**(CpDry/R_d))

#-- 飽和水蒸気圧 e_sat を計算
#-- e_sat = 6.112 exp[{17.67(T-273.15)}/(T-29.65)] [hPa]
         a   = (17.67*(yzTemp-273.15))/(yzTemp-29.65)
         e_sat = 6.112*(a.exp)
#-- 飽和混合比 q_sat を計算
#-- q_sat = (R_d/R_v)(e_sat/p)
         q_sat = (R_d/R_v)*(e_sat/(p-e_sat))

#-- 飽和温度 temp_sat を計算
#-- temp_sat = [1/(T-55)-ln(q_v/q_sat)/2840]^{-1}+55 [k]
         temp_sat = (1.0/(1.0/(yzTemp-55.0)-((yzH2Og/q_sat).log)/2840.0)) + 55.0

#-- 相当温位 theta_e を計算
#-- theta_e = T(Pref/p)**{0.2854(1-0.28 q_v)} *exp[q_v(1+0.81q_v)((3376/T_sat)-2.54)]
         b  = (yzH2Og*(1.0+(0.81*yzH2Og))*((3376.0/temp_sat)-2.54))
         theta_e = yzTemp*(P_0/p)**(0.2854*(1-(0.28*yzH2Og)))*(b.exp)

         theta_e_gphys = theta_e

theta_e_gphys.data.set_att('long_name', "Equivalent potential temerature")
theta_e_gphys.units = 'K'
GPhys::NetCDF_IO.write(ncFile, theta_e_gphys.rename('EPT')) 
# netCDFファイルにライト
#rename で名前を変える
ncFile.close