# -*- 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')
gphysTempAll  = GPhys::NetCDF_IO.open('thermal-moist_TempAll.nc', 'TempAll')
gphysH2OgAll = GPhys::NetCDF_IO.open('thermal-moist_H2O-gAll.nc', 'H2O-gAll')
gphysExnerAll = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll.nc', 'ExnerAll')

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

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

PressBasis = 96500.0e0                # [Pa]
CpDry    = 1039.64234361457397
GasRDry  = 296.785690317661306
MolWtDry = 0.2801348000000000032e-1   # [kg mol-1]
MolWtWet = 0.180152800000000016e-1    # [kg mol-1]


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

gphysLogSVapPress = AMPA / gphysTempAll + AMPB + AMPC * gphysTempAll.log + AMPD * gphysTempAll + AMPE * ( gphysTempAll**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 )  # 確認

p exp( gphysLogSVapPress )  # 確認


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

gphysPressure = gphysExnerAll**( CpDry/GasRDry ) * PressBasis
#gphysHumidity = gphysH2O-gAll / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / PressBasis ) )
gphysHumidity = gphysH2OgAll / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / gphysPressure ) )
p gphysHumidity.val # 確認


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

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

# 出力らへん
outfile = NetCDF.create('thermal-moist_RevHumid.nc')
GPhys::NetCDF_IO.write( outfile, ( gphysHumidity ).rename('RevHumid') )

# outfile を閉じる
outfile.close

exit


# 以下は, 並列計算で rank000000.nc と rank000001.nc が出力されたとき用

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

#gphys = GPhys::NetCDF_IO.open('T.jan.nc', 'T')
gphys1 = GPhys::NetCDF_IO.open('thermal-moist_TempAll_rank000000.nc', 'TempAll')
gphys2 = GPhys::NetCDF_IO.open('thermal-moist_H2O-gAll_rank000000.nc', 'H2O-gAll')
gphys3 = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll_rank000000.nc', 'ExnerAll')

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

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

PressBasis = 96500.0e0                # [Pa]
CpDry    = 1039.64234361457397
GasRDry  = 296.785690317661306
MolWtDry = 0.2801348000000000032e-1   # [kg mol-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 = gphys2 / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / gphysPressure ) )
p gphysHumidity.val # 確認


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

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

# 出力らへん
#outfile = NetCDF.create('tmp.nc')
outfile = NetCDF.create('thermal-moist_RevHumid_rank000000.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('RevHumid') )
#GPhys::NetCDF_IO.write( outfile, gphys )

# outfile を閉じる
outfile.close

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

#gphys = GPhys::NetCDF_IO.open('T.jan.nc', 'T')
gphys1 = GPhys::NetCDF_IO.open('thermal-moist_TempAll_rank000001.nc', 'TempAll')
gphys2 = GPhys::NetCDF_IO.open('thermal-moist_H2O-gAll_rank000001.nc', 'H2O-gAll')
gphys3 = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll_rank000001.nc', 'ExnerAll')

#-- 温度からを飽和蒸気圧を計算 (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 = gphys2 / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / gphysPressure ) )
p gphysHumidity.val # 確認


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

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

# 出力らへん
#outfile = NetCDF.create('tmp.nc')
outfile = NetCDF.create('thermal-moist_RevHumid_rank000001.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('RevHumid') )
#GPhys::NetCDF_IO.write( outfile, gphys )

# outfile を閉じる
outfile.close