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

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

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

#-- 定数

Grav = 9.8
#CpDry = 1039.6423436145739742642
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 水蒸気の気体定数
GasRDry  = 288.17502069917014750899 # deepconv のRd

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

#-- 幅設定

DelX = 500
DelZ =250
DelTime = 150
NT = (21600/150).to_i # どれだけ時間ステップか(NetCDFの軸を生成するために)
NX = 256
NZ = 120

#-- NArray.を用意

zTemp = NArray.sfloat(NX)
p_sat = NArray.sfloat(NX)

#x_fNl = NArray.sfloat(NX)
#x_fAl = NArray.sfloat(NX)

#-- NArray オブジェクト
#-- indgen(0,500)<- 0 から始めて 500 ずつ値を入れる
x_a = VArray.new(NArray.sfloat(NX).indgen(0,500), 
         {"long_name"=>"x-coodinate", "units"=>"m"}, # ハッシュ(ruby)
         "x")


#---- GPhys の軸のオブジェクトに格上げ
x = Axis.new.set_pos(x_a)

z_a = VArray.new(NArray.sfloat(NZ).indgen(0,250), 
         {"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,150),
         {"long_name"=>"Time", "units"=>"sec"}, # ハッシュ(ruby)
          "t" )

#---- 軸のオブジェクトに格上げ
t = Axis.new.set_pos(t_a)

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

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

#-- netCDF づくり
ncFile = NetCDF.create("EPT.nc")
GPhys::NetCDF_IO.write_grid(ncFile, Grid.new(x,z,t))

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

gphysPTemp = GPhys::IO.open('./thermal_PTempAll.nc','PTempAll')
gphysExner = GPhys::IO.open('./thermal_ExnerAll.nc','ExnerAll')
gphysH2Og = GPhys::IO.open('./thermal_H2O-gAll.nc','H2O-gAll')

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

### y
gphysPTemp = gphysPTemp.mean( 'y' )
gphysExner = gphysExner.mean( 'y' )
gphysH2Og = gphysH2Og.mean( 'y' )

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 )
      zPTemp.units = '1'

### NArray オブジェクトにする
      zPTemp = zPTemp.val
      zExner = zExner.val
      zH2Og = zH2Og.val

#-- 温度を計算

      zTemp = zPTemp*zExner

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

#-- 飽和水蒸気圧 e_sat を計算
#-- e_sat = 6.112 exp[{17.67(T-273.15)}/(T-29.65)] [hPa]
         a   = (17.67*(zTemp-273.15))/(zTemp-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)

#-- 飽和温度 temp_sat を計算
#-- temp_sat = [1/(T-55)-ln(q_v/q_sat)/2840]^{-1}+55 [k]
      temp_sat = (1.0/(1.0/(zTemp-55.0)-(log(zH2Og/q_sat))/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  = (zH2Og*(1.0+(0.81*zH2Og))*((3376.0/temp_sat)-2.54))
      theta_e = zTemp*(P_0/p)**(0.2854*(1-(0.28*zH2Og)))*exp(b)

#-- 結果を最終変数に代入
      k = ((z-125)/DelZ).to_i
      n = (t/DelTime).to_i
      theta_e_gphys[0..NX-1,k,n] = theta_e[0..NX-1]
   end
end

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