#! /usr/bin/ruby
# -*- coding: utf-8 -*-
#
# 2013/11/22 石渡正樹
#
# = make_geopotential.rb
# * ジオポテンシャル高度のデータを作成する  未完成!!!!
#
# == USAGE
#   % make_geopotential.rb --time_start=731 --time_end=1096
#
# == 計算方法
# σ座標系における静水圧の式は
#
#      xyz_Height(:,:,1) =                          &
#        & xy_SurfHeight                            &
#        & + GasRDry / Grav * xyz_VirTempWork(:,:,1) * ( 1.0_DP - z_Sigma(1) )
#      do k = 2, kmax
#        xyz_Height(:,:,k) = xyz_Height(:,:,k-1)         &
#          & + GasRDry / Grav * xyr_VirTempWork(:,:,k-1) &
#          &   * r_DelSigma(k-1) / r_Sigma(k-1)
#      end do


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

## 定数設定
GasRUniv = 8.314

MolWtDry = 28.964e-3

GasRDry          = GasRUniv / MolWtDry



EpsV  = MolWtWet / MolWtDry




## オプション解析
parser = GetoptLong.new
parser.set_options(
                   ###    global option   ###
                   ['--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

## データ読み込み, 加工
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')
sigm = GPhys::IO.open('QVap.nc', 'sigm')
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)

sig_weight = GPhys::IO.open('QVap.nc', 'sig_weight')

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

# xyz_VirTemp = xyz_Temp * ( 1.0_DP + ((( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QVap) )


VirTemp = Temp * 


Press_va = VArray.new( NArray.sfloat(nlon, nlat, nsig, ntime),
           {"long_name"=>"pressure", "units"=>"Pa"},
                   "Press" )

for k in 0..nsig-1
  Press_va[true,true,k,true] = na_sig[k] * Ps.val[true,true,true]
end


lon_a = VArray.new( lon.val,
                    {"long_name"=>"longitude", "units"=>"degrees_east"},
                    "lon" )
lon_new = Axis.new.set_pos(lon_a)
 
lat_a = VArray.new( lat.val,
                    {"long_name"=>"latitude","units"=>"degrees_north"},
                    "lat" )
lat_new = Axis.new.set_pos(lat_a)

sig_a = VArray.new( sig.val,
                    {"long_name"=>"sigma at layer midpoints","units"=>"1"},
                    "sig" )
sig_new = Axis.new.set_pos(sig_a)

time_a = VArray.new( time.val,
                    {"long_name"=>"time","units"=>"day"},
                    "time" )
time_new = Axis.new.set_pos(time_a)


Press = GPhys.new( Grid.new(lon_new,lat_new,sig_new,time_new), Press_va)


##### 相対湿度の計算


QVapSat_va = VArray.new( NArray.sfloat(nlon, nlat, nsig, ntime),
           {"long_name"=>"saturation specific humidity", "units"=>"1"},
                   "QVapSat" )

QVapSat_na = EpsV*Es0*NMath::exp(LatentHeat/GasRWet*(1.0/273.0 - 1.0/Temp.val)) / Press.val

QVapSat_va[true,true,true,true] = QVapSat_na[true,true,true,true]

QVapSat =  GPhys.new( Grid.new(lon_new,lat_new,sig_new,time_new), QVapSat_va)


RelHumid_na = QVap.val / QVapSat_na


RelHumid_va = VArray.new( RelHumid_na,
           {"long_name"=>"relative humidity", "units"=>"1"},
                   "RelHumid" )


RelHumid = GPhys.new( Grid.new(lon_new,lat_new,sig_new,time_new), RelHumid_va)


# 書き出し
outfile_qvapsat = NetCDF.create('QVapSat.nc')
outfile_press = NetCDF.create('Press.nc')
outfile_rh = NetCDF.create('RelHumid.nc')

GPhys::NetCDF_IO.write( outfile_qvapsat, QVapSat )
GPhys::NetCDF_IO.write( outfile_press, Press )
GPhys::NetCDF_IO.write( outfile_rh, RelHumid )

outfile_qvapsat.close
outfile_press.close
outfile_rh.close