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

#
# = 概要
#
# * 地表面での水蒸気フラックスを計算するときに用いる
#
#     \frac{\nu_v}{\nu_d} * \frac{e_{sfc}}{p_{sfc}}
#
#   の地表面での水蒸気混合比の term を計算する. 
#   計算結果は thermal-moist_SfcTerm_rank\d\d\d\d\d\d.nc に出力される
# 
# = めも
#
# * 新規作成: 2012/10/02
# * 最終更新: 2012/10/05
#

#-- 必要な netCDF のデータを取得 --#

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]
TempSfc    = 302.0e0                  # 地表面温度       [K]
PressSfc   = 965.0e2                  # 地表面圧力       [Pa]
CpDry    = 1039.64234361457397
GasRDry  = 296.785690317661306
MolWtDry = 0.2801348000000000032e-1   # [kg mol-1]
MolWtWet = 0.180152800000000016e-1    # [kg mol-1]


#-- 地表面温度についての飽和蒸気圧を計算 (AMP の式を使う) し, 結果を配列に入れる --#

LogSVapPress = AMPA / TempSfc + AMPB + AMPC * log( TempSfc ) + AMPD * TempSfc + AMPE * ( TempSfc**2) - log( 1.0e1 )
p LogSVapPress # 確認

SVapPress = exp( LogSVapPress )
p SVapPress  # 確認

# 定数を配列に入れる
gphysSVapPress = ( SVapPress * gphysExnerAll / gphysExnerAll ).cut( 'z'=>150 )
p gphysSVapPress.val # 確認


#-- 最下層のエクスナー関数と標準圧力などから地表面圧力を計算 --#

#gphysPressure = gphys3**( CpDry/GasRDry ) * PressBasis
gphysPressureSfc = ( gphysExnerAll.cut( 'z'=>150 ) )**( CpDry/GasRDry ) * PressBasis
p gphysPressureSfc.val

#-- 地表面での水蒸気混合比を計算 --#

#gphysHumidity = gphys2 / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / PressBasis ) )
gphysH2OgAllSfc =  ( MolWtWet / MolWtDry ) * ( gphysSVapPress / gphysPressureSfc ) 
p gphysH2OgAllSfc.val # 確認

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

# 変数の rename など
#gphysHumidity.data.rename!("RevHumid")
gphysH2OgAllSfc.data.rename!("H2O-gAll_Sfc")
#gphysHumidity.data.set_att('long_name', "rerative humidity")
gphysH2OgAllSfc.data.set_att('long_name', "H2O-g Mixing Ratio at Surface")

# 出力らへん
outfile = NetCDF.create('thermal-moist_H2O-gAll_Sfc.nc')
GPhys::NetCDF_IO.write( outfile, ( gphysH2OgAllSfc ).rename('H2O-gAll_Sfc') )

# outfile を閉じる
outfile.close

exit

#######################################################
##### thermal-moist_H2O-gSfc_rank000001.nc を作成 #####
#######################################################

#-- 必要な netCDF のデータを取得 --#

#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')
gphysExnerAll = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll_rank000001.nc', 'ExnerAll')

#-- 地表面温度についての飽和蒸気圧を計算 (AMP の式を使う) し, 結果を配列に入れる --#

LogSVapPress = AMPA / TempSfc + AMPB + AMPC * log( TempSfc ) + AMPD * TempSfc + AMPE * ( TempSfc**2) - log( 1.0e1 )
p LogSVapPress # 確認

SVapPress = exp( LogSVapPress )
p SVapPress  # 確認

# 定数を配列に入れる
gphysSVapPress = ( SVapPress * gphysExnerAll / gphysExnerAll ).cut( 'z'=>150 )
p gphysSVapPress.val # 確認


#-- 最下層のエクスナー関数と標準圧力などから地表面圧力を計算 --#

#gphysPressure = gphys3**( CpDry/GasRDry ) * PressBasis
gphysPressureSfc = ( gphysExnerAll.cut( 'z'=>150 ) )**( CpDry/GasRDry ) * PressBasis
p gphysPressureSfc.val

#-- 地表面での水蒸気混合比を計算 --#

#gphysHumidity = gphys2 / ( ( MolWtWet / MolWtDry ) * ( gphysSVapPress / PressBasis ) )
gphysH2OgAllSfc =  ( MolWtWet / MolWtDry ) * ( gphysSVapPress / gphysPressureSfc ) 
p gphysH2OgAllSfc.val # 確認

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

# 変数の rename など
#gphysHumidity.data.rename!("RevHumid")
gphysH2OgAllSfc.data.rename!("H2O-gAll_Sfc")
#gphysHumidity.data.set_att('long_name', "rerative humidity")
gphysH2OgAllSfc.data.set_att('long_name', "H2O-g Mixing Ratio at Surface")

# 出力らへん
#outfile = NetCDF.create('thermal-moist_RevHumid_rank000001.nc')
outfile = NetCDF.create('thermal-moist_H2O-gAll_Sfc_rank000001.nc')
#GPhys::NetCDF_IO.write( outfile, ( gphysHumidity ).rename('RevHumid') )
GPhys::NetCDF_IO.write( outfile, ( gphysH2OgAllSfc ).rename('H2O-gAll_Sfc') )
#GPhys::NetCDF_IO.write( outfile, gphys )

# outfile を閉じる
outfile.close