#! /usr/bin/ruby
# -*- coding: utf-8 -*-
#
# 2007/09/28 石渡正樹
# 2013/05/26 石渡正樹
#
# = make_RH.rb
# * 相対湿度のデータを作成する
#
# == USAGE
#   % make_RH.rb --saturation_scheme=nha --time_start=731 --time_end=1095
#   % make_RH.rb --saturation_scheme=nha --time_start=1 --time_end=3
# /GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/make_RH_v2.rb --saturation_scheme=tetens --time_start=1 --time_end=2
# /GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/make_RH_v2.rb --saturation_scheme=nha 

require 'getoptlong'
require "numru/ggraph"
include NumRu

## 定数設定
GasRUniv = 8.314
LatentHeat = 2.5e6
# latent heat for NHA1992 [J mol-1]
LatHeat = 43655.0
MolWtWet = 18.01528e-3
MolWtDry = 28.964e-3
Es0 =611.0

GasRWet          = GasRUniv / MolWtWet
EpsV  = MolWtWet / MolWtDry

P0 = 1.4e+11

## オプション解析
parser = GetoptLong.new
parser.set_options(
                   ###    global option   ###
                   ['--saturation_scheme',                GetoptLong::REQUIRED_ARGUMENT],
                   ['--lon',                GetoptLong::REQUIRED_ARGUMENT],
                   ['--time_start',         GetoptLong::REQUIRED_ARGUMENT],
                   ['--time_end',           GetoptLong::REQUIRED_ARGUMENT],
                   ['--range',              GetoptLong::REQUIRED_ARGUMENT],
                   ['--wsn',                GetoptLong::REQUIRED_ARGUMENT]
                   )
parser.each_option do |name, arg|
  eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'"  # strage option value to $OPT_val
end

time_start = ($OPT_time_start||1).to_f
time_end = ($OPT_time_end||9999).to_f
saturation_scheme = ($OPT_saturation_scheme||'nha')


## データ読み込み, 加工
QVap = GPhys::IO.open('QVap.nc', 'QVap').cut('time'=>time_start..time_end)
Temp = GPhys::IO.open('Temp.nc', 'Temp').cut('time'=>time_start..time_end)
Ps = GPhys::IO.open('Ps.nc', 'Ps').cut('time'=>time_start..time_end)

sig = GPhys::IO.open('QVap.nc', 'sig')
lat = GPhys::IO.open('QVap.nc', 'lat')
lon = GPhys::IO.open('QVap.nc', 'lon')
#time = GPhys::IO.open('QVap.nc', 'time').cut('time'=>time_start..time_end)

na_sig = sig.val
na_lon = lon.val
na_lat = lat.val
#na_time = time.val

nsig = na_sig.size
nlon = na_lon.size
nlat = na_lat.size
#ntime = na_time.size


dir = './'

ofile = NetCDF.create(dir + 'RelHumid.nc')

# 時間ループ
GPhys::NetCDF_IO.each_along_dims_write([Ps,Temp,QVap], ofile, 'time') { 
  |ps,temp,qvap|  
        
#  time = qvap.axis("time")    

  rh = temp.copy

  # 気圧の計算
  na_ps = ps.val
  na_press = temp.copy
  na_press = na_ps.reshape(nlon, nlat ,1, 1) * na_sig.reshape(1, 1 ,nsig, 1)

  # 飽和比湿の計算
  if saturation_scheme == 'nha' then # (NHA1992) の式
    na_qvsat = EpsV* (P0/na_press) * NMath.exp(-LatHeat/(GasRUniv*temp.val))
  elsif saturation_scheme == 'tetens' then # Tetens の式
    na_qvsat = EpsV*Es0*NMath.exp(LatentHeat/GasRWet*(1.0/273.0 - 1.0/temp.val))/na_press
  else
    p 'SATURATION SCHEME IS INVALID!!' 
    exit
  end

  # 相対湿度の計算
  rh.val = na_qvsat
  rh = qvap/rh

  rh.units = '1/1'
  rh.long_name = 'relative humidity'
  rh.name = 'RelHumid'

  [rh]
}
ofile.close