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

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

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

#-- 定数

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

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

#-- 幅設定

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

#-- 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_a = VArray.new(NArray.sfloat(NY).indgen(0,DelY), 
         {"long_name"=>"y-coodinate", "units"=>"m"}, # ハッシュ(ruby)
         "y")
y = Axis.new.set_pos(y_a)

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)

#-- 物理場オブジェクト
theta_es_a = VArray.new(NArray.sfloat(NX,NY,NZ,NT+1),
         {"long_name"=>"saturation equivalent potential temperature", "units"=>"(1)"},
         "theta_es")

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

#-- netCDF づくり
ncFile = NetCDF.create("satEPT.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 )

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

### NArray オブジェクトにする
         yzPTemp = yzPTemp.val
         yzExner = yzExner.val
         yzH2Og = yzH2Og.val

#-- 温度を計算

         yzTemp = yzPTemp*yzExner

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

#-- 飽和水蒸気圧 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*exp(a)
#-- 飽和混合比 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(yzH2Og/q_sat))/2840.0)) + 55.0

#-- 相当温位 theta_es を計算
#-- theta_es = T(Pref/p)**{0.2854(1-0.28 q_v)} *exp[q_v(1+0.81q_v)((3376/T_sat)-2.54)]
#-- These are only replaced yzH2Og in cal_EPT.rb to q_sat.  
         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)))*exp(b)

#-- 結果を最終変数に代入
         l = ((y-(DelY/2).to_i)/DelY).to_i
         k = ((z-(DelZ/2).to_i)/DelZ).to_i
         n = (t/DelTime).to_i
         theta_es_gphys[0..NX-1,l,k,n] = theta_es[0..NX-1]
      end
   end
end

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('theta_es')) 
# netCDFファイルにライト
#rename で名前を変える
ncFile.close