#! /usr/bin/ruby
# -*- coding: euc-jp -*-
#
# 2013-5-12 石渡正樹
# OLR-Ts 関係 : 1 次元放射対流平衡モデル, GCM

require "getoptlong"        # for option_parse
require "numru/ggraph"
include NumRu


# 2 次元以上のデータを 1 次元配列に落とす (GGraph.scatter 対策)
def flatten_gphys(gp)
  na_dat = gp.val.flatten
  n = na_dat.size
  axis = dummy_axis(n)
  data = VArray.new( na_dat,
                     {"long_name"=>gp.data.long_name, 
                       "units"=>gp.data.units.to_s},
                     gp.data.name )
  gp = GPhys.new( Grid.new(axis), data )
end

# ダミーの gphys を作成 (GGraph.scatter での図の範囲指定のため)
def dummy_gphys(na)
  n = na.size
  axis = dummy_axis(n)
  data = VArray.new( na,
                     {"long_name"=>"dummy", 
                       "units"=>"1"},
                     "dummy" )
  gp = GPhys.new( Grid.new(axis), data )
end

def dummy_axis(n)
  axis_a = VArray.new( NArray.sfloat(n).indgen,
                       {"long_name"=>"dummy","units"=>"1"},
                       "dummy" )
  axis = Axis.new.set_pos(axis_a)
end

# 領域平均. ただし lon のみ cut したものを使っても構わない
def ave_area(gp, na_lat_weight)
  gp = gp.mean('lon') * na_lat_weight
  gp = 0.5 * gp.sum('lat')
  gp.val = [gp.val]
end

##########################################################################

## オプション解析
parser = GetoptLong.new
parser.set_options(
  ###    global option   ###
  ['--area',                      GetoptLong::REQUIRED_ARGUMENT],
  ['--time_start',                GetoptLong::REQUIRED_ARGUMENT],
  ['--time_end',                  GetoptLong::REQUIRED_ARGUMENT],
  ['--onedim',                      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

#データ格納ディレクトリ
DATA_HOME_ONED = "/GFD_Dennou_Work3/momoko/OneDimEquilibrium"
GCMDATA_HOME = "/GFD_Dennou_Work10/momoko/SyncRotEarthRad-2"
GCMDIR_HEADER = "SR_CLT0_Omega"
GCMDIR_FOOTER = "_T42L26"


#GCMDIR_HEADER = "SR_CLT1500_Omega"
#GCMDIR_FOOTER = "_S1800_T42L26"


#ディレクトリリスト
oned_dirs = ['Earth_ConstStrat_L32_rh100', 'Earth_ConstStrat_L32_rh80', 'Earth_ConstStrat_L32_rh60', 'Earth_ConstStrat_L32_rh40']

time_start = ($OPT_time_start||1).to_f
time_end = ($OPT_time_end||9999).to_f

area = ($OPT_area||'180').to_s

onedim=($OPT_onedim||'on').to_s

#omegas = [0.0, 0.1, 0.5, 1.0]
omegas = [0.5]
nomega = omegas.size

#RH_SCATTER_MAX = 1.0
#RH_SCATTER_MIN = 0.4
SCATTER_MAX = 3
SCATTER_MIN = 0

RPlanet = 6.371 * 10 ** 6  # [m]

TITLE = 'OLR'
YNAME = 'OLR'
#YUNIT = '10^15 W'

longname_rh='relative humidity'
name_rh= 'RH'
units_rh= '1'

P_0 = 1.4 * 1.0e11          # NHA92 で用いる基準気圧 [Pa]
Grav = 9.8                  # 重力加速度 [m s-2]
GasRUniv = 8.314            # 普遍気体定数 [J K-1 mol-1]
MolWtWet = 18.01528e-3      # 凝結成分の平均分子量 [kg mol-1]
MolWtDry = 28.964e-3        # 乾燥成分の平均分子量 [kg mol-1]

XMIN = 200
XMAX = 400
YMIN = 120
YMAX = 420

#XMIN = 250
#XMAX = 320
#YMIN = 100
#YMAX = 330


DCL.gropn($OPT_wsn||4)

DCL.sgpset('lfull', true)     # 全画面表示
DCL.sgpset('lclip', true)     # 図の外側は切り取り
DCL.sgpset('lcntl', true) 
DCL.uzfact(0.7)
GGraph.set_fig( 'itr'=> 1, 'viewport'=>[0.2,0.8,0.3,0.6] )

#
# 1 次元データの 1 次元グラフ
#

if onedim == 'on' then
  for i in 0..oned_dirs.size-1
    ## データ切出し, Ts 軸データの作成
    olr = GPhys::IO.open("#{DATA_HOME_ONED}/#{oned_dirs[i]}/OLRB.nc", 'OLRB').cut('time'=>0.0,'lon'=>0.0)
    units_olr = 'W m|-2"'
    olr_a = VArray.new( olr.val, 
                      {"long name"=>"outgoing longwave","units"=>units_olr},
                      "OLR")

    surftemp = GPhys::IO.open("#{DATA_HOME_ONED}/#{oned_dirs[i]}/sst_1dim.nc", 
                    'SurfTemp').cut('time'=>0.0,'lon'=>0.0)
    surftemp_a = VArray.new(surftemp.val,
              {"long name"=>"Surface Temperature","units"=>"K"},"SurfTemp")
    surftemp_axis = Axis.new.set_pos(surftemp_a)

    olr_temp_axis = GPhys.new(Grid.new(surftemp_axis), olr_a)

    ## 作図
    if i == 0 then

# 範囲指定のためのダミーデータ
na_xrange = NArray.to_na([XMIN, XMAX])
na_yrange = NArray.to_na([YMIN, YMAX])
gp_xrange = dummy_gphys(na_xrange)
gp_yrange = dummy_gphys(na_yrange)
gp_xrange.data.long_name = surftemp.data.long_name
gp_yrange.data.long_name = olr.data.long_name
gp_xrange.data.units = surftemp.data.units
gp_yrange.data.units = olr_temp_axis.data.units

# 図の範囲を指定する方法が分からないのでダミーの点を打つ
GGraph.scatter( gp_xrange, gp_yrange, true, 'type'=>1)

#      GGraph.line( olr_temp_axis, false, "min"=>YMIN, "max"=>YMAX, "index"=>855, "annotate"=>false, "exchange"=>false )
      GGraph.line( olr_temp_axis, false, "min"=>YMIN, "max"=>YMAX, "index"=>5, "annotate"=>false, "exchange"=>false )
    else
#      index = 855 - i*150 # 色つき
      index = 5
      p index
      GGraph.line( olr_temp_axis, false, "min"=>YMIN, "max"=>YMAX, "index"=>index, "annotate"=>false, "exchange"=>false )
    end
  end # end of oned data loop
end

#
# GCM のデータの散布図
#

hash_gp_Tg = Hash.new
hash_gp_OLR = Hash.new

## 各 omega で回す.
for i in 0..nomega-1

  omega = omegas[i]

  dir = GCMDIR_HEADER + omega.to_s + GCMDIR_FOOTER

  ### 読み込み
  # OLR and basic value
  var = 'OLRA'
  file = var + '.nc'
  path = File.join(GCMDATA_HOME, dir, file)
  gp_OLR = GPhys::IO.open(path, var)

  sig = GPhys::IO.open(path, 'sig')
  sig_weight = GPhys::IO.open(path, 'sig_weight')
  na_lon = GPhys::IO.open(path, 'lon').val
  na_lat = GPhys::IO.open(path, 'lat').val
  na_lat_weight = GPhys::IO.open(path, 'lat_weight').val
  nlon = na_lon.size
  nlat = na_lat.size
  nsig = sig.val.size

  # SurfTemp
  var = 'SurfTemp'
  file = var + '.nc'
  path = File.join(GCMDATA_HOME, dir, file)
  gp_Tg = GPhys::IO.open(path, var)

  if (area == '360') then # 全球領域
    hash_gp_OLR[omega] = gp_OLR.cut('time'=>time_start..time_end).mean('time')
    hash_gp_Tg[omega]  = gp_Tg.cut('time'=>time_start..time_end).mean('time')
  elsif (area == '180') then  # 直下点を中心とする 180度×180度領域
    hash_gp_OLR[omega] = gp_OLR.cut('time'=>time_start..time_end,'lon'=>0..180, 'lat'=>-90..90).mean('time')
    hash_gp_Tg[omega]  = gp_Tg.cut('time'=>time_start..time_end,'lon'=>0..180, 'lat'=>-90..90).mean('time')
  elsif (area == '120') then # 直下点を中心とする 120度×120度領域
    hash_gp_OLR[omega] = gp_OLR.cut('time'=>time_start..time_end,'lon'=>30..150, 'lat'=>-60..60).mean('time')
    hash_gp_Tg[omega]  = gp_Tg.cut('time'=>time_start..time_end,'lon'=>30..150, 'lat'=>-60..60).mean('time')
  elsif (area == '90') then # 直下点を中心とする 90度×90度領域
    hash_gp_OLR[omega] = gp_OLR.cut('time'=>time_start..time_end,'lon'=>45..135, 'lat'=>-45..45).mean('time')
    hash_gp_Tg[omega]  = gp_Tg.cut('time'=>time_start..time_end,'lon'=>45..135, 'lat'=>-45..45).mean('time')
  else
    p 'area is not set properly: area=', area
    exit
  end

end # omega のループ終わり



opt = {
  'max'=>YMAX, 'min'=>YMIN,  
  'annotate'=>false,
  'title'=>TITLE,
  'index'=>1
}


opt['type'] = 1  # 4
opt['min'] = SCATTER_MIN
opt['max'] = SCATTER_MAX

for i in 0..omegas.size-1
  omega = omegas[i]
  gp_Tg = flatten_gphys(hash_gp_Tg[omega])
  gp_OLR = flatten_gphys(hash_gp_OLR[omega])
  gp_color = gp_Tg.copy
  gp_color.val = gp_color.val.fill(i)
  if !(onedim == 'on') && i == 0 then 
    GGraph.color_scatter( gp_Tg, gp_OLR, gp_color, true, opt)
  else
    GGraph.color_scatter( gp_Tg, gp_OLR, gp_color, 
                      false, opt)
  end
end

DCL.grcls


#-85.76059: 250K
#-80.26878, -74.74454, -69.21297, -63.67863, -58.14296: 300K
#-52.60653, -47.06964, -41.53246, -35.99508, -30.45755: 350K
#-24.91993, -19.38223, -13.84448, -8.306703, -2.768903: 400K
# 2.768903,  8.306703,  13.84448,  19.38223,  24.91993: 450K
# 30.45755,  35.99508,  41.53246,  47.06964,  52.60653: 500K
# 58.14296,  63.67863,  69.21297,  74.74454,  80.26878: 550K
# 85.76059: 560K