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

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("satEPT.nc")

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

gphysPTemp = GPhys::IO.open('./thermal1_PTempAll.nc','PTempAll')
gphysExner = GPhys::IO.open('./thermal1_ExnerAll.nc','ExnerAll')

### 時間でカット
   yzPTemp = gphysPTemp
   yzExner = gphysExner

#-- 温度を計算
         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)-(log(1))/2840.0)) + 55.0

#-- 飽和相当温位 theta_es を計算(H2Ogに飽和混合比を代入)
#-- theta_es = T(Pref/p)**{0.2854(1-0.28 q_v)} *exp[q_v(1+0.81q_v)((3376/T_sat)-2.54)]
         b  = (q_sat*(1.0+(0.81*q_sat))*((3376.0/temp_sat)-2.54))
         theta_es = yzTemp*(P_0/p)**(0.2854*(1-(0.28*q_sat)))*(b.exp)

         theta_es_gphys = theta_es

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