#!/usr/bin/env ruby

# ----------------------------------------------
# local load path

$local_path = '/home/yukiko/work/ape/yukiko/lib' 
$: << $local_path

# ----------------------------------------------
# 必要なライブラリ, モジュールの読み込み

load "#{$local_path}/ape-view.rb"
load "#{$local_path}/lib-ape-view.rb"
load "#{$local_path}/lib-ape-base.rb"

# --------------------------------------------------------------------


END{


  ["mradlAc", "con", "mradlAa", "mradlAb", "mradlAd" ].each{ |expid|
    $expid = expid
    $cum = "adj"
    gr_comp_mknetcdf
    $cum = "kuo"
    gr_comp_mknetcdf

  }

  
}


# ----------------------------------------------

def set_selection_point(method)

  if $expid == "con" && $cum == "adj" 
    if $filter == "ad"
      threshold = 0.0002 
    elsif $filter == "k"
      threshold = 0.0001  
    elsif $filter == "g"
      threshold = 0.0001  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAa" && $cum == "adj" 
    if $filter == "ad"
      threshold = 0.0002  
    elsif $filter == "k"
      threshold = 0.0001
    elsif $filter == "g"
      threshold = 0.0001
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAb" && $cum == "adj" 
    if $filter == "ad"
      threshold = 0.0002  
    elsif $filter == "k"
      threshold = 0.0001  
    elsif $filter == "g"
      threshold = 0.0001  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAc" && $cum == "adj" 
    if $filter == "ad"
      threshold = 0.0003  
    elsif $filter == "k"
      threshold = 0.0002  
    elsif $filter == "g"
      threshold = 0.0002  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAd" && $cum == "adj" 
    if $filter == "ad"
      threshold = 0.0003  
    elsif $filter == "k"
      threshold = 0.0002  
    elsif $filter == "g"
      threshold = 0.0002  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "con" && $cum == "kuo" 
    if $filter == "ad"
      threshold = 0.0002  
    elsif $filter == "k"
      threshold = 0.00005  
    elsif $filter == "g"
      threshold = 0.00005  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAa" && $cum == "kuo" 
    if $filter == "ad"
      threshold = 0.0002  
    elsif $filter == "k"
      threshold = 0.00005  
    elsif $filter == "g"
      threshold = 0.00005  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAb" && $cum == "kuo" 
    if $filter == "ad"
      threshold = 0.0002  
    elsif $filter == "k"
      threshold = 0.00005  
    elsif $filter == "g"
      threshold = 0.00005  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAc" && $cum == "kuo" 
    if $filter == "ad"
      threshold = 0.0001  
    elsif $filter == "k"
      threshold = 0.00005  
    elsif $filter == "g"
      threshold = 0.00005  
    elsif $filter == "n"
      threshold = 0.0012
    end
  elsif $expid == "mradlAd" && $cum == "kuo" 
    if $filter == "ad"
      threshold = 0.0002  
    elsif $filter == "k"
      threshold = 0.0002  
    elsif $filter == "g"
      threshold = 0.00005  
    elsif $filter == "n"
      threshold = 0.0012
    end
  end
  return threshold
end

# ----------------------------------------------

def gr_comp_mknetcdf(method="threshold")

  filt_array = ["ad","k","g"]
  filt_hash =  { 
    "ad" => "advect", 
    "k"  => "kelvin",
    "g"  => "grav"
  }
  
  filt_array.each{ |filter|
    
    
    $filter = filter
    threshold = set_selection_point(method)
    lost_axes = "comp threshold = #{threshold}"  

    $outfile = "AGCM5-#{$cum}-#{$expid}_#{filt_hash[filter]}filtreal-composite.nc"

    dir_path = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/spctfilt/"
    @data =
      Ape.new("#{dir_path}AGCM5-#{$cum}-#{$expid}_#{filt_hash[filter]}filt.nc")

    @data_real =
      Ape.new("#{dir_path}AGCM5-#{$cum}-#{$expid}_nofilt.nc")

    
    #  tr_tppn
    rain = @data.gphys_open("tr_tppn_#{filter}filt")
    rain = rain.add_lost_axes(rain.get_att("lost_axis"))
    rain = rain[true,-401..-1]
    rain = rain.
      set_lost_axes(rain.lost_axes.to_s ).
      add_lost_axes("refer to #{filter} spct filtering data").
      add_lost_axes(lost_axes)
  
    # comp_point composite
    gphys = @data_real.gphys_open("tr_t")
    gphys = gphys[true,true,-401..-1]
    comp = gphys.composite(rain, threshold, true)
    lost =
      "#{gphys.get_att("lost_axis")}\nrefer to #{filt_hash[filter]} spct filtering data\n#{lost_axes}"
    comp = comp.      
      set_att("ape_name",
              "composite_point" ).
      set_att("units", "1" ).
      rename("comp_point").
      set_att("lost_axis",lost)
  
    lost_diff = lost + "\n(diff) from (mean) zonal"

    # T composite
    gphys = @data_real.gphys_open("tr_t")
    gphys = gphys[true,true,-401..-1]
    t_comp = gphys.composite(rain, threshold)
    t_comp = t_comp.
      set_att("ape_name", "t_(u,_-sigdot)_composite").
      set_att("units", "K, (m s-1, s-1)" )
    t_comp = (t_comp - t_comp.mean(0)).
      set_att("lost_axis",lost_diff). 
      rename("comp_t_#{filter}filtreal")
  
    # Q composite
    gphys = @data_real.gphys_open("tr_q")
    gphys = gphys[true,true,-401..-1]
    q_comp = gphys.composite(rain, threshold)
    q_comp = q_comp.
      set_att("ape_name", "q_(u,_-sigdot)_composite").
      set_att("units", "1, (m s-1, s-1)" )
    q_comp = (q_comp - q_comp.mean(0)).
      set_att("lost_axis",lost_diff).
      rename("comp_q_#{filter}filtreal")
    
    # tconv composite
    gphys = @data_real.gphys_open("tr_tconv")
    gphys = gphys[true,true,-401..-1]
    tconv_comp = gphys.composite(rain, threshold)
    tconv_comp = tconv_comp.
      set_att("ape_name", "qcnd_(u,_-sigdot)_composite").
      set_att("units", "K s-1, (m s-1, s-1)" )
    tconv_comp = tconv_comp.
      set_att("lost_axis",lost).
      rename("comp_tconv_#{filter}filtreal")
    
    # U composite
    gphys = @data_real.gphys_open("tr_u")
    gphys = gphys[true,true,-401..-1]
    u_comp = gphys.composite(rain, threshold)
    u_comp = (u_comp - u_comp.mean(0)).
      set_att("ape_name", "u_(u,_-sigdot)_composite").
      set_att("units", "m s-1, (m s-1, s-1)" ).
      set_att("lost_axis",lost_diff).
      rename("comp_u_#{filter}filtreal")
          
    # V composite
    gphys = @data_real.gphys_open("tr_v")
    gphys = gphys[true,true,-401..-1]
    v_comp = gphys.composite(rain, threshold)
    v_comp = (v_comp - v_comp.mean(0)).
      set_att("ape_name", "v_(u,_-sigdot)_composite").
      set_att("units", "m s-1, (m s-1, s-1)" ).
      set_att("lost_axis",lost_diff).
      rename("comp_v_#{filter}filtreal")

    # OMG composite
    gphys = @data_real.gphys_open("tr_om")
    gphys = gphys[true,true,-401..-1]
    om_comp = gphys.composite(rain, threshold)    
    om_comp = ( om_comp - om_comp.mean(0) ).
      set_att("lost_axis",lost_diff).
      rename("comp_om_#{filter}filtreal")

    # Gphys オブジェクト生成  & nc ファイル生成  ----
    
    # 出力ファイル
    ncfile = NetCDF.create("#{$outfile}")
    
    # composite point
    GPhys::NetCDF_IO.write( ncfile, comp )
    
    # composite した変数
    GPhys::NetCDF_IO.write( ncfile, t_comp )
    GPhys::NetCDF_IO.write( ncfile, q_comp )
    GPhys::NetCDF_IO.write( ncfile, tconv_comp )
    GPhys::NetCDF_IO.write( ncfile, u_comp )
    GPhys::NetCDF_IO.write( ncfile, v_comp )
    GPhys::NetCDF_IO.write( ncfile, om_comp )
  
    ncfile.close  
  
  }

end

# ---------------------------------------------------------------

class GPhys

  def composite(rain, threshold, flag=false )

    if threshold.class == Float
      composite_threshold(rain, threshold, flag)
    else
      composite_online(rain, threshold, flag )
    end

  end

  def composite_threshold(rain, threshold, flag=false )

    # gphys : コンポジットをとる変数 (lon,sigmap,time)
    # rain  : コンポジットを評価する rain データ (lon,time)
    # flag=true  : コンポジットの参照点 (lon,time) を返す
    # flag=false : gphys のコンポジットの結果 (lon,sigmap) を返す
    # threshold  : rain 閾値
 
    # 閾値
    # threshold = 0.0005
    
    # 軸の値を保存
    grid = @grid.copy
    grid_lon = grid.axis(0)
    grid_sigmap = grid.axis(1)
    grid_time = grid.axis(2)
    
    # 軸のサイズを保存
    grid_lon_size    = grid.axis(0).pos.val.size
    grid_sigmap_size = grid.axis(1).pos.val.size
    grid_time_size   = grid.axis(2).pos.val.size

    # 作業領域初期化 
    rain_cal = 
      NArray.sfloat( rain.coord(0).val.size + 6 ).fill!(0.0)
    comp_num = 0 ; count = 0
    comp_point =  NArray.sfloat(grid_lon_size,grid_time_size).fill!(0.0)
    
    # データ
    data = @data.copy.val
    
    # コンポジット値初期化
    comp = 
      NArray.sfloat(grid_lon_size , grid_sigmap_size ).fill!(0.0)

    grid_time_size.times{ |time|
    
      # コンポジット格子値計算用 rain データ作成
      rain_cal[3..-4] = rain.val[true,time]
      
      # 重ね合わせ
      (rain_cal.size-6).times{ |comp_num|
	
	# 周囲よりも値が大きければその点を中心経度に移動
	if rain_cal[comp_num+3] > threshold then
	  num = comp_num + 3
	  
	  if  rain_cal[num] > rain_cal[num - 3] && 
	      rain_cal[num] > rain_cal[num - 2] && 
	      rain_cal[num] > rain_cal[num - 1] && 
	      rain_cal[num] > rain_cal[num + 1] && 
	      rain_cal[num] > rain_cal[num + 2] && 
	      rain_cal[num] > rain_cal[num + 3]    then 
          
	    comp_point[comp_num,time] = 1.0
	    
	    print '(lgrid,tgrid) = (', comp_num, ',', time, ')', "\n"
	    puts rain_cal[num]
	    
	    n =  grid_lon_size/2 - comp_num
	    
	    tmp = NArray.sfloat(grid_lon_size,grid_sigmap_size).fill!(0.0)
	    nm = grid_lon_size


	    if n == 0 then
	      
	      tmp = data[true,true,time]
	      
	    elsif n > 0 then
	      
	      tmp[0..(n-1),true]  = data[(nm-n)..(nm-1),true,time] 
	      tmp[n..(nm-1),true]    = data[0..(nm-1-n),true,time]
	      
	    else 
	      
	      tmp[(nm+n)..(nm-1),true]  = data[0..(-1-n),true,time] 
	      tmp[0..(nm-1+n),true]    = data[(-n)..(nm-1),true,time]

	    end  

	    comp = comp + tmp  
	    count = count + 1 
	    
	  end	  

	end
	
      }

    }  

    # 平均
    comp = comp / count

    # Gphys オブジェクト生成  ----
    
    # composite point
    grid = Grid.new(grid_lon, grid_time)
    comp_point = VArray.new(comp_point).rename("comp_point")
    comp_point = GPhys.new(grid, comp_point) 
    
    # composite した変数
    grid = Grid.new(grid_lon, grid_sigmap)
    comp = VArray.new(comp).rename("#{@data.copy.name}_comp")
    comp = GPhys.new(grid, comp)
    
    if flag 
      return comp_point
    else
      return comp
    end

  end
   
end