#!/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|
  ["mradlAc" ].each{ |expid|
    $expid = expid
    $cum = "adj"
    gr_comp_mknetcdf
#    $cum = "kuo"
#    gr_comp_mknetcdf

  }
    
}


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

def set_selection_point(method)
  
  if $expid == "mradlAc" && $cum == "adj" 
    if $filter == "km"
      threshold = 7.0
    elsif $filter == "kf"
      threshold = 2.5
    end
  end
  return threshold
end

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

def gr_comp_mknetcdf(method="threshold")

#  filt_array = ["ad","k","g","n"]
#  filt_array = ["n"]
  filt_array = ["kf","km"]
  filt_hash =  { 
    "km"  => "moist_kelvin",
    "kf"  => "free_kelvin"
  }
  
  filt_array.each{ |filter|
    
    $filter = filter
    threshold = set_selection_point(method)
    lost_axes = "comp div(u_z12) threshold = #{threshold}"  

    $outfile = "AGCM5-#{$cum}-#{$expid}_#{filter}filt-u_z12-composite.nc"

    dir_path = "/home/yukiko/work/aqua3d/yukiko/sh/ruby/mradlstdexp/mlshtr/work-AGCM5-spctfilt/test/"

    @data =
      Ape.new("#{dir_path}AGCM5-#{$cum}-#{$expid}_#{filter}filt.nc")
    
    #  tr_tppn
    if filter == "n"
      rain = @data.gphys_open("tr_u")
    else
      rain = @data.gphys_open("tr_u_#{filter}filt")
    end

    rain = rain.add_lost_axes(rain.get_att("lost_axis"))
    rain = rain[true,11,-401..-1]
    rain = rain.
      set_lost_axes(rain.lost_axes.to_s).
      add_lost_axes("#{filter} spct filtering").
      add_lost_axes(lost_axes)

    rain = rain.sabun
    
    # comp_point composite
    if filter == "n"
      gphys = @data.gphys_open("tr_t")
    else
      gphys = @data.gphys_open("tr_t_#{filter}filt")
    end
    gphys = gphys.
      add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma"))
    gphys = gphys[true,true,-401..-1]
    lost =
      "#{gphys.lost_axes.to_s.sub("y=","lat=")}\n#{filt_hash[filter]} spct filtering\n#{lost_axes}"
    comp = gphys.composite(rain, threshold, true)
    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
    if filter == "n"
      gphys = @data.gphys_open("tr_t")
    else
      gphys = @data.gphys_open("tr_t_#{filter}filt")
    end
    gphys = gphys.
      add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma"))
    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}filt")
  
    # Q composite
    if filter == "n"
      gphys = @data.gphys_open("tr_q")
    else
      gphys = @data.gphys_open("tr_q_#{filter}filt")
    end
    gphys = gphys.
      add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma"))
    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}filt")
  
    # tr_tconv composite
    if filter == "n"
      gphys = @data.gphys_open("tr_tconv")
    else
      gphys = @data.gphys_open("tr_tconv_#{filter}filt")
    end
    gphys = gphys.
      add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma"))
    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}filt")

    # U composite
    if filter == "n"
      gphys = @data.gphys_open("tr_u")
    else
      gphys = @data.gphys_open("tr_u_#{filter}filt")
    end
    gphys = gphys.
      add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma"))
    gphys = gphys[true,true,-401..-1]
    u_comp = gphys.composite(rain, threshold)
    u_comp = u_comp.
      set_att("ape_name", "u_(u,_-sigdot)_composite").
      set_att("units", "m s-1, (m s-1, s-1)" )
    u_comp = (u_comp - u_comp.mean(0)).
      set_att("lost_axis",lost_diff).
      rename("comp_u_#{filter}filt")

    # V composite
    if filter == "n"
      gphys = @data.gphys_open("tr_v")
    else
      gphys = @data.gphys_open("tr_v_#{filter}filt")
    end
    gphys = gphys.
      add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma"))
    gphys = gphys[true,true,-401..-1]
    v_comp = gphys.composite(rain, threshold)
    v_comp = v_comp.
      set_att("ape_name", "v_(u,_-sigdot)_composite").
      set_att("units", "m s-1, (m s-1, s-1)" )
    v_comp = (v_comp - v_comp.mean(0)).
      set_att("lost_axis",lost_diff).
      rename("comp_v_#{filter}filt")

    # OMG composite
    if filter == "n"
      gphys = @data.gphys_open("tr_om")
    else
      gphys = @data.gphys_open("tr_om_#{filter}filt")
    end
    gphys = gphys.
      add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma"))
    gphys = gphys[true,true,-401..-1]
    om_comp = gphys.composite(rain, threshold)
    om_comp = om_comp.
      set_att("ape_name", "sigdot_(u,_-sigdot)_composite").
      set_att("units", "s-1, (m s-1, s-1)" )
    om_comp = (om_comp - om_comp.mean(0)).
      set_att("lost_axis",lost_diff).
      rename("comp_om_#{filter}filt")
    
    # 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)
    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

            tmp[n..-1,true] = data[0..-1-n,true,time]
            tmp[0..n-1,true] = data[-n..-1,true,time]

	    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

  def sabun
    x = self.dup.copy
    y = self.dup.copy
    x[1..-1,true] = self[0..-2,true]
    x[0,true] = self[-1,true]
    y[-1,true] = self[0,true]
    y[0..-2,true] = self[1..-1,true]
    return (y-x)
  end

  def cshift_x(n)
    y = self.dup
    y[n..-1,true] = self[0..-1-n,true]
    y[0..n-1,true] = self[-n..-1,true]
    return y
  end

  def cshift_y(n)
    y = self.dup
    y[true,n..-1] = self[true,0..-1-n]
    y[true,0..n-1] = self[true,-n..-1]
    return y
  end

end