#!/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_horizontal_mknetcdf
    $cum = "kuo"
    gr_comp_horizontal_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_horizontal_mknetcdf(method="threshold")

  $horizontal = true
  
  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-horizontal.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-horizontal.nc")

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

    # 出力ファイル
    ncfile = NetCDF.create("#{$outfile}")
    
    lost =
      "refer to #{filt_hash[filter]} spct filtering data\n#{lost_axes}"
    lost_diff = lost + "\n(diff) from (mean) zonal"


    #  tr_tppn
    rain = @data.gphys_open("tr_tppn_#{filter}filt").cut(true,0,true)
    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 )
    
    # comp_point composite
    gphys = @data_real.gphys_open("tr_tppn") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    comp = gphys.composite(rain, threshold, true)
    lost =
      "#{rain.lost_axes.to_s}\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)
    GPhys::NetCDF_IO.write( ncfile, comp )
    
    # tr_tppn composite
    gphys = @data_real.gphys_open("tr_tppn") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    tppn_comp = gphys.composite(rain, threshold)
    tppn_comp = tppn_comp.
      set_att("ape_name", "rain_(u,_v)_composite").
      set_att("units", "kg m-2 s-1, (m s-1, m s-1)" )
    tppn_comp = tppn_comp.
      set_att("lost_axis",lost). 
      rename("comp_tppn_#{filter}filtreal")

      GPhys::NetCDF_IO.write( ncfile, tppn_comp )
  
    # tr_mslp composite
    gphys = @data_real.gphys_open("tr_mslp") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    mslp_comp = gphys.composite(rain, threshold)
    mslp_comp = mslp_comp.
      set_att("ape_name", "ps_(u,_v)_composite").
      set_att("units", "Pa, (m s-1, m s-1)" )
    mslp_comp = (mslp_comp - mslp_comp.mean(0)).
      set_att("lost_axis",lost_diff). 
      rename("comp_mslp_#{filter}filtreal")
    
    GPhys::NetCDF_IO.write( ncfile, mslp_comp )


    # tr_u_z12 composite
    gphys = @data_real.gphys_open("tr_u_z12") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    u250s_comp = gphys.composite(rain, threshold)
    u250s_comp = u250s_comp.
      set_att("ape_name", "u_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    u250s_comp = (u250s_comp - u250s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_u_z12_#{filter}filtreal")
    
    GPhys::NetCDF_IO.write( ncfile, u250s_comp )

    # tr_u_z8 composite
    gphys = @data_real.gphys_open("tr_u_z8") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    u570s_comp = gphys.composite(rain, threshold)
    u570s_comp = u570s_comp.
      set_att("ape_name", "u_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    u570s_comp = (u570s_comp - u570s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_u_z8_#{filter}filtreal")

      GPhys::NetCDF_IO.write( ncfile, u570s_comp )

    # tr_u_z5 composite
    gphys = @data_real.gphys_open("tr_u_z5") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    u850s_comp = gphys.composite(rain, threshold)
    u850s_comp = u850s_comp.
      set_att("ape_name", "u_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    u850s_comp = (u850s_comp - u850s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_u_z5_#{filter}filtreal")
    
    GPhys::NetCDF_IO.write( ncfile, u850s_comp )

    # tr_u_z3 composite
    gphys = @data_real.gphys_open("tr_u_z3") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    u950s_comp = gphys.composite(rain, threshold)
    u950s_comp = u950s_comp.
      set_att("ape_name", "u_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    u950s_comp = (u950s_comp - u950s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_u_z3_#{filter}filtreal")
    
    GPhys::NetCDF_IO.write( ncfile, u950s_comp )

    # tr_v_z12 composite
    gphys = @data_real.gphys_open("tr_v_z12") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    v250s_comp = gphys.composite(rain, threshold)
    v250s_comp = v250s_comp.
      set_att("ape_name", "v_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    v250s_comp = (v250s_comp - v250s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_v_z12_#{filter}filtreal")
    
    GPhys::NetCDF_IO.write( ncfile, v250s_comp )

    # tr_v_z8 composite
    gphys = @data_real.gphys_open("tr_v_z8") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    v570s_comp = gphys.composite(rain, threshold)
    v570s_comp = v570s_comp.
      set_att("ape_name", "v_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    v570s_comp = (v570s_comp - v570s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_v_z8_#{filter}filtreal")

    GPhys::NetCDF_IO.write( ncfile, v570s_comp )

    # tr_v_z5 composite
    gphys = @data_real.gphys_open("tr_v_z5") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true) 
    v850s_comp = gphys.composite(rain, threshold)
    v850s_comp = v850s_comp.
      set_att("ape_name", "v_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    v850s_comp = (v850s_comp - v850s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_v_z5_#{filter}filtreal")

    GPhys::NetCDF_IO.write( ncfile, v850s_comp )

    # tr_v_z3 composite
    gphys = @data_real.gphys_open("tr_v_z3") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    v950s_comp = gphys.composite(rain, threshold)
    v950s_comp = v950s_comp.
      set_att("ape_name", "v_(u,_v)_composite").
      set_att("units", "m s-1, (m s-1, m s-1)" )
    v950s_comp = (v950s_comp - v950s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_v_z3_#{filter}filtreal")
    
    GPhys::NetCDF_IO.write( ncfile, v950s_comp )

    # tr_gph_z12 composite
    gphys = @data_real.gphys_open("tr_gph_z12") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    phi250s_comp = gphys.composite(rain, threshold)
    phi250s_comp = phi250s_comp.
      set_att("ape_name", "gph_(u,_v)_composite").
      set_att("units", "m, (m s-1, m s-1)" )
    phi250s_comp = (phi250s_comp - phi250s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_gph_z12_#{filter}filtreal")

    GPhys::NetCDF_IO.write( ncfile, phi250s_comp )

    # tr_gph_z5 composite
    gphys = @data_real.gphys_open("tr_gph_z5") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    phi850s_comp = gphys.composite(rain, threshold)
    phi850s_comp = phi850s_comp.
      set_att("ape_name", "gph_(u,_v)_composite").
      set_att("units", "m, (m s-1, m s-1)" )
    phi850s_comp = (phi850s_comp - phi850s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_gph_z5_#{filter}filtreal")

    GPhys::NetCDF_IO.write( ncfile, phi850s_comp )

    # tr_t570s composite
    gphys = @data_real.gphys_open("tr_t_z8") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    t570s_comp = gphys.composite(rain, threshold)
    t570s_comp = t570s_comp.
      set_att("ape_name", "t_(u,_v)_composite").
      set_att("units", "K, (m s-1, m s-1)" )
    t570s_comp = (t570s_comp - t570s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_t_z8_#{filter}filtreal")
    
    GPhys::NetCDF_IO.write( ncfile, t570s_comp )

    # tr_q_z5 composite
    gphys = @data_real.gphys_open("tr_q_z5") 
    gphys = gphys[true,true,-401..-1].
      cut(true,-30..30,true)
    q850s_comp = gphys.composite(rain, threshold)
    q850s_comp = q850s_comp.
      set_att("ape_name", "q_(u,_v)_composite").
      set_att("units", "1, (m s-1, m s-1)" )
    q850s_comp = (q850s_comp - q850s_comp.mean(0)).
      set_att("lost_axis","#{gphys.get_att("lost_axis")}\n#{lost_diff}"). 
      rename("comp_q_z5_#{filter}filtreal")

    GPhys::NetCDF_IO.write( ncfile, q850s_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