# -*- coding: euc-jp -*-
require "numru/gphys"
include NumRu

#
# = 概要
#
# * 飽和水蒸気圧を計算して, その結果を用いて相対湿度を計算する
#   計算結果は thermal-moist_RevHumid_rank\d\d\d\d\d\d.nc に出力される
#

##### thermal-moist_ReHumid.nc を作成 #####

#gphys = GPhys::NetCDF_IO.open('T.jan.nc', 'T')
gphys1 = GPhys::NetCDF_IO.open('thermal1_PTempAll.nc', 'PTempAll')
gphys2 = GPhys::NetCDF_IO.open('thermal1_H2O-gAll.nc', 'H2O-gAll')
gphys3 = GPhys::NetCDF_IO.open('thermal1_ExnerAll.nc', 'ExnerAll')
gphys1 = gphys1*gphys3 # 温度を計算

#-- 相対湿度の計算に必要な定数を設定 --#

AMPA = -2313.0338e0
AMPB = -164.03307e0
AMPC =  38.053682e0
AMPD = -0.13844344e0
AMPE =  7.4465367e-5

PressBasis = 101000.0e0                # [Pa]
CpDry    = 1004
GasRDry  = 287.0
#MolWtDry = 0.2801348000000000032e-1   # [kg mol-1]
MolWtDry = 0.2896864e-1
MolWtWet = 0.180152800000000016e-1    # [kg mol-1]

#-- 温度からを飽和蒸気圧を計算 (AMP の式を使う) --#

gphysLogSVapPress = AMPA / gphys1 + AMPB + AMPC * gphys1.log + AMPD * gphys1 + AMPE * ( gphys1**2) - log( 1.0e1 )
#p gphys1.val             # 確認
#p gphysLogSVapPress.val  # 確認

gphysSVapPress = gphysLogSVapPress.exp
p gphysSVapPress.val  # 確認

gphysLogSVapPress = AMPA / 373.15 + AMPB + AMPC * log( 373.15 ) + AMPD * 373.15 + AMPE * ( 373.15** 2) - log( 1.0e1 )  # 確認
#gphysLogSVapPress2 = AMPA / 298.58 + AMPB + AMPC * log( 298.58 ) + AMPD * 298.58 + AMPE * ( 298.58** 2) - log( 1.0e1 ) # 確認

p exp( gphysLogSVapPress )  # 確認
#p exp( gphysLogSVapPress2 ) # 確認


#-- さっき求めた飽和水上気圧と水蒸気混合比などから相対湿度を計算する --#

gphysPressure = gphys3**( CpDry/GasRDry ) * PressBasis
#gphysHumidity = gphys2 / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / PressBasis ) )
gphysHumidity = 100.0*gphys2 / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / gphysPressure ) )
p gphysHumidity.val # 確認


#-- NetCDF ファイルへの書き込み --#

# 変数の rename など
gphysHumidity.data.rename!("Hum")
gphysHumidity.data.set_att('long_name', "relative humidity")

# 出力らへん
#outfile = NetCDF.create('tmp.nc')
outfile = NetCDF.create('thermal1_Hum.nc')
#GPhys::NetCDF_IO.write( outfile, gphys.cut('level'=>1000..250).mean(0) )
#GPhys::NetCDF_IO.write( outfile, gphys.cut('level'=>1000..250).mean(0,1).rename('T00') )
GPhys::NetCDF_IO.write( outfile, ( gphysHumidity ).rename('Hum') )
#GPhys::NetCDF_IO.write( outfile, gphys )

# outfile を閉じる
outfile.close