#!/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{

  $rezol = "T39L48_eml"
#  gr_comp_mknetcdf("line")
#  gr_comp_mknetcdf
  gr_comp_horizontal_mknetcdf("line")
  gr_comp_horizontal_mknetcdf
  (1..7).each{ |num|
#    gr_comp_mknetcdf("line#{num}")
    gr_comp_horizontal_mknetcdf("line#{num}")
  }

  $rezol = "T39L48_ias"
#  gr_comp_mknetcdf("line")
#  gr_comp_mknetcdf  
  gr_comp_horizontal_mknetcdf("line")
  gr_comp_horizontal_mknetcdf
  (1..3).each{ |num|
#    gr_comp_mknetcdf("line#{num}")
    gr_comp_horizontal_mknetcdf("line#{num}")
  }

  ["ksc","kuo","mca","non"].each{ |item|

    $rezol = "T39L48_#{item}"
#    gr_comp_mknetcdf
#    gr_comp_mknetcdf("line")
  gr_comp_horizontal_mknetcdf("line")
  gr_comp_horizontal_mknetcdf
  }

  ["eml","non"].each{ |item|

    $rezol = "T79L48_#{item}"
#    gr_comp_mknetcdf
#    gr_comp_mknetcdf("line")
  gr_comp_horizontal_mknetcdf("line")
  gr_comp_horizontal_mknetcdf

    $rezol = "T39L96_#{item}"
#    gr_comp_mknetcdf
#    gr_comp_mknetcdf("line")
  gr_comp_horizontal_mknetcdf("line")
  gr_comp_horizontal_mknetcdf

    $rezol = "T39L24_#{item}"
#    gr_comp_mknetcdf
#    gr_comp_mknetcdf("line")
  gr_comp_horizontal_mknetcdf("line")
  gr_comp_horizontal_mknetcdf

#    $rezol = "T159L48_#{item}"
##    gr_comp_mknetcdf
#    gr_comp_mknetcdf("line")

  }


#  ["T39L96_non","T159L48_eml","T159L48_non"].each{ |rezol|
  ["T39L96_non"].each{ |rezol|
    $rezol = rezol
#    gr_comp_mknetcdf("line1")
  gr_comp_horizontal_mknetcdf("line1")
  }


    $rezol = "T159L48_eml"
##    gr_comp_mknetcdf
#    gr_comp_mknetcdf("line")
#    gr_comp_mknetcdf("line1")
#  gr_comp_horizontal_mknetcdf("line")
#  gr_comp_horizontal_mknetcdf("line1")
#  gr_comp_horizontal_mknetcdf

    $rezol = "T159L48_non"
#  gr_comp_horizontal_mknetcdf("line")
#  gr_comp_horizontal_mknetcdf("line1")
#  gr_comp_horizontal_mknetcdf
  
}

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

def gr_comp_horizontal_mknetcdf(method="threshold")

  $horizontal = true

  $gr2ncfile_path = "/home/yukiko/work/ape/yukiko/data/#{$rezol}_expID01/"
  $outfile = "/home/yukiko/work/ape/yukiko/data/composite/#{$rezol}_expID01_composite_horizontal_#{method}.nc"

  threshold = set_selection_point(method)
  if threshold.class == Float
    lost_axes_add = "comp threshold = #{threshold}"  
  else
    lost_axes_add = "comp (lon,time)=(#{threshold["lonmin"]},#{threshold["tmin"]})-(#{threshold["lonmax"]},#{threshold["tmax"]})"
  end
  lost_diff = "#{lost_axes_add}\n(diff) from (mean) zonal"
  Dir.foreach($gr2ncfile_path) { |item|
    if item =~ /expID01_horizontal_00/
      @data = Ape.new("#{$gr2ncfile_path}/#{item}") 
    end
  }


  if $rezol == "T159L48_eml"
    rain = __gr_tr_tppn[true,-801..-401]
  else
    rain = __gr_tr_tppn[true,-401..-1]
  end
  lost_axes = rain.lost_axes.to_s.sub("y=","lat=")
  lost_axes << "\n#{lost_axes_add}"
  rain = rain.set_lost_axes(lost_axes).lon_lotate

  # comp_point composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_u250s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_u250s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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_axes)
  
  # tppn composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_tppn")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_tppn")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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_axes_add}"). 
    rename("comp_tppn")
  
  # mslp composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_mslp")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_mslp")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_mslp")

  # u250s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_u250s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_u250s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_u250s")

  p "u250s #{gphys.get_att("lost_axis")}"

  # u570s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_u570s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_u570s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_u570s")

  p "u570s #{gphys.get_att("lost_axis")}"

  # u850s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_u850s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_u850s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_u850s")

  p "u850s #{gphys.get_att("lost_axis")}"

  # u950s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_u950s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_u950s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_u950s")

  p "u950s #{gphys.get_att("lost_axis")}"

  # v250s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_v250s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_v250s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_v250s")

  p "v250s #{gphys.get_att("lost_axis")}"

  # v570s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_v570s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_v570s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_v570s")

  p "v570s #{gphys.get_att("lost_axis")}"

  # v850s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_v850s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_v850s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_v850s")

  p "v850s #{gphys.get_att("lost_axis")}"

  # v950s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_v950s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_v950s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_v950s")

  p "v950s #{gphys.get_att("lost_axis")}"

  # phi250s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_phi250s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_phi250s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  phi250s_comp = gphys.composite(rain, threshold)
  phi250s_comp = phi250s_comp.
    set_att("ape_name", "PHI_(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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_phi250s")

  p "phi250s #{gphys.get_att("lost_axis")}"

  # phi850s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_phi850s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_phi850s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  phi850s_comp = gphys.composite(rain, threshold)
  phi850s_comp = phi850s_comp.
    set_att("ape_name", "PHI_(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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_phi850s")

  p "phi850s #{gphys.get_att("lost_axis")}"

  # t570s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_t570s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_t570s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_t570s")

  p "t570s #{gphys.get_att("lost_axis")}"

  # q850s composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("tr_q850s")[true,true,-801..-401].
      cut(true,-30..30,true).lon_lotate
  else
    gphys = __gr_tr_lat0("tr_q850s")[true,true,-401..-1].
      cut(true,-30..30,true).lon_lotate
  end
  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}"). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_q850s")

  p "q850s #{gphys.get_att("lost_axis")}"

  # Gphys オブジェクト生成  & nc ファイル生成  ----
  
  # 出力ファイル
  ncfile = NetCDF.create("#{$filepath_out}/#{$outfile}")

  # composite point
  GPhys::NetCDF_IO.write( ncfile, comp )
  
  # composite した変数
  GPhys::NetCDF_IO.write( ncfile, tppn_comp )
  GPhys::NetCDF_IO.write( ncfile, mslp_comp )
  GPhys::NetCDF_IO.write( ncfile, u250s_comp )
  GPhys::NetCDF_IO.write( ncfile, u570s_comp )
  GPhys::NetCDF_IO.write( ncfile, u850s_comp )
  GPhys::NetCDF_IO.write( ncfile, u950s_comp )
  GPhys::NetCDF_IO.write( ncfile, v250s_comp )
  GPhys::NetCDF_IO.write( ncfile, v570s_comp )
  GPhys::NetCDF_IO.write( ncfile, v850s_comp )
  GPhys::NetCDF_IO.write( ncfile, v950s_comp )
  GPhys::NetCDF_IO.write( ncfile, phi250s_comp )
  GPhys::NetCDF_IO.write( ncfile, phi850s_comp )
  GPhys::NetCDF_IO.write( ncfile, t570s_comp )
  GPhys::NetCDF_IO.write( ncfile, q850s_comp )
  
  ncfile.close  
  
  $horizontal = false

end

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

def gr_comp_mknetcdf(method="threshold")

  $gr2ncfile_path = "/home/yukiko/work/ape/yukiko/data/#{$rezol}_expID01/"
  $outfile = "/home/yukiko/work/ape/yukiko/data/composite/#{$rezol}_expID01_composite_#{method}.nc"
  threshold = set_selection_point(method)
  if threshold.class == Float
    lost_axes_add = "\ncomp threshold = #{threshold}"  
  else
    lost_axes_add = "\ncomp (lon,time)=(#{threshold["lonmin"]},#{threshold["tmin"]})-(#{threshold["lonmax"]},#{threshold["tmax"]})"
  end
  Dir.foreach($gr2ncfile_path) { |item|
    if item =~ /expID01_00/
      @data = Ape.new("#{$gr2ncfile_path}/#{item}") 
    end
  }


  if $rezol == "T159L48_eml"
    rain = __gr_tr_tppn[true,-801..-401]
  else
    rain = __gr_tr_tppn[true,-401..-1]
  end
  lost_axes = rain.lost_axes.to_s.sub("y=","lat=")
  lost_axes << lost_axes_add
  rain = rain.set_lost_axes(lost_axes).lon_lotate

  # comp_point composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("T")[true,true,-801..-401].lon_lotate
  else
    gphys = __gr_tr_lat0("T")[true,true,-401..-1].lon_lotate
  end
  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_axes)
  
  # T composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("T")[true,true,-801..-401].lon_lotate
  else
    gphys = __gr_tr_lat0("T")[true,true,-401..-1].lon_lotate
  end
  t_comp = gphys.composite(rain, threshold)
  t_comp = t_comp.
    set_att("ape_name", "T_(U,_-OMG)_composite").
    set_att("units", "K, (m s-1, Pa s-1)" )
  t_comp = (t_comp - t_comp.mean(0)).
    set_att("lost_axis",lost_axes). 
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_t")
  
  # Tv composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("TV")[true,true,-801..-401].lon_lotate
  else
    gphys = __gr_tr_lat0("TV")[true,true,-401..-1].lon_lotate
  end
  tv_comp = gphys.composite(rain, threshold)
  tv_comp = tv_comp.
    set_att("ape_name", "Tv_(U,_-OMG)_composite").
    set_att("units", "K, (m s-1, Pa s-1)" )
  tv_comp = (tv_comp - tv_comp.mean(0)).
    set_att("lost_axis",lost_axes).
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_tv")
  
  # Q composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("Q")[true,true,-801..-401].lon_lotate
  else
    gphys = __gr_tr_lat0("Q")[true,true,-401..-1].lon_lotate
  end
  q_comp = gphys.composite(rain, threshold)
  q_comp = q_comp.
    set_att("ape_name", "Q_(U,_-OMG)_composite").
    set_att("units", "1, (m s-1, Pa s-1)" )
  q_comp = (q_comp - q_comp.mean(0)).
    set_att("lost_axis",lost_axes).
    add_lost_axes("(diff) from (mean) zonal").
    rename("comp_q")
  
  # tconv composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("pf_t_conv")[true,true,-801..-401].lon_lotate
  else
    gphys = __gr_tr_lat0("pf_t_conv")[true,true,-401..-1].lon_lotate
  end
  tconv_comp = gphys.composite(rain, threshold)
  tconv_comp = tconv_comp.
    set_att("ape_name", "DTCOND_(U,_-OMG)_composite").
    set_att("units", "K s-1, (m s-1, Pa s-1)" )
  tconv_comp = tconv_comp.
    set_att("lost_axis",lost_axes).
    rename("comp_tconv")
  
  # U composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("U")[true,true,-801..-401].lon_lotate
  else
    gphys = __gr_tr_lat0("U")[true,true,-401..-1].lon_lotate
  end
  u_comp = gphys.composite(rain, threshold)
  u_comp = (u_comp - u_comp.mean(0)).
    set_att("ape_name", "U_(U,_-OMG)_composite").
    set_att("units", "m s-1, (m s-1, Pa s-1)" ).
    set_att("lost_axis",lost_axes).
    rename("comp_u")
  
  # OMG composite
  if $rezol == "T159L48_eml"
    gphys = __gr_tr_lat0("OMG")[true,true,-801..-401].lon_lotate
  else
    gphys = __gr_tr_lat0("OMG")[true,true,-401..-1].lon_lotate
  end
  om_comp = gphys.composite(rain, threshold)    
  om_comp = ( om_comp - om_comp.mean(0) ).
    set_att("lost_axis",lost_axes).
    rename("comp_om")  

  # 軸の long_name 直し (pressure level => sigma)
  z = t_comp.grid_copy.axis(1).pos.
    set_att("long_name","sigma").
    set_att("units","1")
  z = Axis.new().set_pos(z)
  grid = t_comp.grid_copy.change_axis(1, z)
  t_comp = GPhys.new( grid, t_comp.data )
  tv_comp = GPhys.new( grid, tv_comp.data )
  q_comp = GPhys.new( grid, q_comp.data )
  tconv_comp = GPhys.new( grid, tconv_comp.data )
  u_comp = GPhys.new( grid, u_comp.data )
  om_comp = GPhys.new( grid, om_comp.data )
  
  
  # Gphys オブジェクト生成  & nc ファイル生成  ----
  
  # 出力ファイル
  ncfile = NetCDF.create("#{$filepath_out}/#{$outfile}")
  
  # composite point
  GPhys::NetCDF_IO.write( ncfile, comp )
  
  # composite した変数
  GPhys::NetCDF_IO.write( ncfile, t_comp )
  GPhys::NetCDF_IO.write( ncfile, tv_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, om_comp )
  
  ncfile.close  
  
end

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

def set_selection_point(method)

  $threshold_array = nil

  if $rezol == "T39L48_eml"
    if method == "line"
      threshold_hash = {
        "tmin" => 45,
        "tmax" => 49,
        "lonmin" => 30, 
        "lonmax" => 100
      }
    elsif method == "line1"
      threshold_hash = {
        "tmin" => 43.5,
        "tmax" => 48.5,
        "lonmin" => -120, 
        "lonmax" => -60
      }
    elsif method == "line2"
      threshold_hash = {
        "tmin" => 75,
        "tmax" => 76.5,
        "lonmin" => -40, 
        "lonmax" => -20
      }
    elsif method == "line3"
      threshold_hash = {
        "tmin" => 63.5,
        "tmax" => 41,
        "lonmin" => -120, 
        "lonmax" => -40
      }
    elsif method == "line4"
      threshold_hash = {
        "tmin" => 74.5,
        "tmax" => 76.5,
        "lonmin" => -30, 
        "lonmax" => 20
      }
    elsif method == "line5"
      threshold_hash = {
        "tmin" => 38,
        "tmax" => 70,
        "lonmin" => -180, 
        "lonmax" => 180
      }
      threshold1_hash = {
        "tmin" => 70,
        "tmax" => 100,
        "lonmin" => -180, 
        "lonmax" => 157.5
      }
      threshold2_hash = {
        "tmin" => 6,
        "tmax" => 38,
        "lonmin" => -180, 
        "lonmax" => 180
      }
      threshold3_hash = {
        "tmin" => 0,
        "tmax" => 6,
        "lonmin" => 112.5, 
        "lonmax" => 180
      }
      $threshold_array = 
        [ threshold1_hash, threshold2_hash, threshold3_hash ]
    elsif method == "line6"
      threshold_hash = {
        "tmin" => 35,
        "tmax" => 55,
        "lonmin" => -180, 
        "lonmax" => 180
      }
      threshold1_hash = {
        "tmin" => 55,
        "tmax" => 75,
        "lonmin" => -180, 
        "lonmax" => 180
      }
      threshold2_hash = {
        "tmin" => 75,
        "tmax" => 95,
        "lonmin" => -180, 
        "lonmax" => 180
      }
      threshold3_hash = {
        "tmin" => 95,
        "tmax" => 100,
        "lonmin" => -180, 
        "lonmax" => -90
      }
      threshold4_hash = {
        "tmin" => 15,
        "tmax" => 35,
        "lonmin" => -180, 
        "lonmax" => 180
      }
      threshold5_hash = {
        "tmin" => 0,
        "tmax" => 15,
        "lonmin" => -90, 
        "lonmax" => 180
      }
      $threshold_array = 
        [ threshold1_hash, threshold2_hash, threshold3_hash,
        threshold4_hash, threshold5_hash ]
    elsif method == "line7"
      threshold_hash = {
        "tmin" => 35,
        "tmax" => 75,
        "lonmin" => -180, 
        "lonmax" => 180
      }
      threshold1_hash = {
        "tmin" => 75,
        "tmax" => 100,
        "lonmin" => -180, 
        "lonmax" => 45
      }
      threshold2_hash = {
        "tmin" => 0,
        "tmax" => 35,
        "lonmin" => -143, 
        "lonmax" => 180
      }
      $threshold_array = 
        [ threshold1_hash, threshold2_hash ]
    end

    threshold = 0.0007  
  elsif $rezol == "T39L48_ias"

    if method == "line"
#      threshold_hash = {
#        "tmin" => 30,
#        "tmax" => 40,
#        "lonmin" => -100, 
#        "lonmax" => -10
#      }
      threshold_hash = {
        "tmin" => 50,
        "tmax" => 55,
        "lonmin" => 10, 
        "lonmax" => 47
      }
    elsif method == "line1"
      threshold_hash = {
        "tmin" => 15,
        "tmax" => 20,
        "lonmin" => -120, 
        "lonmax" => -75
      }
    elsif method == "line2"
      threshold_hash = {
        "tmin" => 64.5,
        "tmax" => 67.5,
        "lonmin" => -150, 
        "lonmax" => -110
      }
    elsif method == "line3"
      threshold_hash = {
        "tmin" => 47,
        "tmax" => 51,
        "lonmin" => -50, 
        "lonmax" => 0
      }
    end

    threshold = 0.0007  
  elsif $rezol == "T39L48_ksc"
    threshold_hash = {
      "tmin" => 48,
      "tmax" => 51,
      "lonmin" => 30, 
      "lonmax" => 75
    }
    threshold = 0.0010  
  elsif $rezol == "T39L48_kuo"
    threshold_hash = {
      "tmin" => 19,
      "tmax" => 24,
      "lonmin" => -150, 
      "lonmax" => -70
    }
    threshold = 0.0010 
  elsif $rezol == "T39L48_mca"

    threshold_hash = {
      "tmin" => 44.5,
      "tmax" => 51.5,
      "lonmin" => 80, 
      "lonmax" => 170
    }

    threshold = 0.0010  
  elsif $rezol == "T39L48_non"
    threshold_hash = {
      "tmin" => 26.5,
      "tmax" => 31.5,
      "lonmin" => -70, 
      "lonmax" => -10
    }
    threshold = 0.0010  
  elsif $rezol == "T39L24_eml"
    threshold_hash = {
      "tmin" => 0,
      "tmax" => 100,
      "lonmin" => -20, 
      "lonmax" => 70
    }
    threshold = 0.0007  
  elsif $rezol == "T39L96_eml"
    threshold_hash = {
      "tmin" => 23.5,
      "tmax" => 28.5,
      "lonmin" => -30, 
      "lonmax" => 30
    }
    threshold = 0.0005  
  elsif $rezol == "T79L48_eml"
    threshold_hash = {
      "tmin" => 79,
      "tmax" => 83,
      "lonmin" => 50, 
      "lonmax" => 150
    }
    threshold = 0.0007  
  elsif $rezol == "T159L48_eml"

    if method == "line"
      threshold_hash = {
#        "tmin" => 73.5,
#        "tmax" => 77.5,
        "tmin" => 73,
        "tmax" => 76.5,
        "lonmin" => 0, 
        "lonmax" => 60
#        "lonmax" => 100
      }
    elsif method == "line1"
      # 2006/01/14 いぜん
      threshold_hash = {
        "tmin" => 74.8,
        "tmax" => 82.56,
        "lonmin" => 0, 
        "lonmax" => 100
      }
    end

    threshold = 0.0007  
  elsif $rezol == "T39L24_non"
    threshold_hash = {
      "tmin" => 26.5,
      "tmax" => 31.5,
      "lonmin" => 70, 
      "lonmax" => 120
    }
    threshold = 0.0010  
  elsif $rezol == "T39L96_non"

    if method == "line"
      threshold_hash = {
        "tmin" => 24,
        "tmax" => 29,
        "lonmin" => 30, 
        "lonmax" => 80
      }
    elsif method == "line1"
      threshold_hash = {
        "tmin" => 60,
        "tmax" => 65.5,
        "lonmin" => 100, 
        "lonmax" => 170
      }
    end
    threshold = 0.0010  
  elsif $rezol == "T79L48_non"
    threshold_hash = {
      "tmin" => 65,
      "tmax" => 69.5,
      "lonmin" => -150, 
      "lonmax" => -100
    }
    threshold = 0.0010  
  elsif $rezol == "T159L48_non"
    if method == "line"
      threshold_hash = {
        "tmin" => 62,
        "tmax" => 71,
        "lonmin" => -50, 
        "lonmax" => 50
      }
    elsif method == "line1"
      # 2006/01/14 いぜん
      threshold_hash = {
        "tmin" => 62.9,
        "tmax" => 72.08,
        "lonmin" => -50, 
        "lonmax" => 50
      }
    end
    threshold = 0.0010  
  end

  if method =~ /line/
    threshold = threshold_hash
  end

  return threshold
end

# --------------------------------------------------------------------
def __gr_tr_tppn
  g = @data.gphys_open("tr_tppn").cut(true,0,true)
  g = g.set_lost_axes(g.lost_axes.to_s.sub("y=","lat=") )
  gphys = __gr2_gphys_make_2d(g,0)
  return gphys
end

def __gr_tr_lat0(var)
  g = @data.gphys_open(var)
  lost_axis = [g.data.get_att("lost_axis") ,
    g.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
  g = g.set_lost_axes(lost_axis)
  
  gphys = __gr2_gphys_make_3d(g). 
  rename("tr_#{var.downcase.gsub("omg","om").gsub("z","phi")}")
  return gphys
end

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

def __gr2_gphys_make_2d(g,dimcut)
  
  grid0size = g.grid_copy.axis(0).pos.val.size
  grid_0    = g.grid_copy.axis(0)
  timesize = 360*4
  grid_time =
    Axis.new().
    set_pos(VArray.new(NArray.sfloat(timesize).indgen!/4+ 0.25).rename("time").
      put_att("units","days"))

  data =  NArray.sfloat(grid0size,timesize).fill!(0.0)
  grid = Grid.new(grid_0, grid_time)

  if $horizontal == true
   hori = "_horizontal"
  else
   hori = ""
  end

  if $rezol =~ /T159L48/

    data[true,0..(3*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0011.nc",
                          g.name).cut(true,dimcut,true).data.val
    data[true,(3*30*4)..(3*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0012.nc",
                          g.name).cut(true,dimcut,true).data.val
    data[true,(3*30*4*2)..(3*30*4*3-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0013.nc",
                          g.name).cut(true,dimcut,true).data.val
    data[true,(3*30*4*3)..(3*30*4*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0014.nc",
                          g.name).cut(true,dimcut,true).data.val
  else

    data[true,0..(6*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0006.nc",
                          g.name).cut(true,dimcut,true).data.val
    data[true,(6*30*4)..(6*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0007.nc",
                          g.name).cut(true,dimcut,true).data.val
  end

  data = VArray.new(data).rename(g.name)
  gphys = GPhys.new(grid, data)
  gphys = gphys.set_att("ape_name",g.get_att("ape_name")). 
                set_att("units",g.get_att("units")). 
                set_lost_axes(g.lost_axes)
#                add_lost_axes(g.get_att("lost_axis"))
  return gphys

end


def __gr2_gphys_make_3d(g)


  if $horizontal == true
   hori = "_horizontal"
  else
   hori = ""
  end

  grid0size = g.grid_copy.axis(0).pos.val.size
  grid_0    = g.grid_copy.axis(0)

  grid1size = g.grid_copy.axis(1).pos.val.size
  grid_1    = g.grid_copy.axis(1)

  timesize = 360*4
  grid_time =
    Axis.new().
    set_pos(VArray.new(NArray.sfloat(timesize).indgen!/4+ 0.25).rename("time").
      put_att("units","days"))

  data =  NArray.sfloat(grid0size,grid1size,timesize).fill!(0.0)
  grid = Grid.new(grid_0, grid_1, grid_time)

  if $rezol =~ /T159L48/

    data[true,true,0..(3*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0011.nc",
                          g.name).data.val
    data[true,true,(3*30*4)..(3*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0012.nc",
                          g.name).data.val
    data[true,true,(3*30*4*2)..(3*30*4*3-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0013.nc",
                          g.name).data.val
    data[true,true,(3*30*4*3)..(3*30*4*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0014.nc",
                          g.name).data.val
  else

    data[true,true,0..(6*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0006.nc",
                          g.name).data.val
    data[true,true,(6*30*4)..(6*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01#{hori}_0007.nc",
                          g.name).data.val
  end

  data = VArray.new(data).rename(g.name)
  gphys = GPhys.new(grid, data)
  gphys = gphys.set_att("ape_name",g.get_att("ape_name")). 
                set_att("units",g.get_att("units")). 
                set_att("lost_axis",g.get_att("lost_axis")). 
                set_lost_axes(g.lost_axes)
#                add_lost_axes(g.get_att("lost_axis"))

  return gphys

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

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

  def composite_online(rain, threshold_hash, flag=false )
 
    # 軸の値を保存
    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
    
    # 作業領域初期化 
    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)

    tmin = threshold_hash["tmin"]
    tmax = threshold_hash["tmax"]
    lonmin = threshold_hash["lonmin"]
    lonmax = threshold_hash["lonmax"]

    grid_lon.pos.val.each{ |comp_num|
      if comp_num > lonmin
	if comp_num < lonmax
          time = 
            tmin*1.0 + 
            (tmax -tmin)*1.0 / ( lonmax -lonmin ) * (comp_num -lonmin)
	  comp,comp_point,count = 
	    composite_calc(time,comp_num,data,comp,comp_point,count)
	end
      end
    }

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

    if $threshold_array.class == Array

    $threshold_array.size.times{ |num| 

      threshold_hash = $threshold_array[num]

      tmin = threshold_hash["tmin"]
      tmax = threshold_hash["tmax"]
      lonmin = threshold_hash["lonmin"]
      lonmax = threshold_hash["lonmax"]

      grid_lon.pos.val.each{ |comp_num|
        if comp_num > lonmin
	  if comp_num < lonmax
            time = 
              tmin*1.0 + 
              (tmax -tmin)*1.0 / ( lonmax -lonmin ) * (comp_num -lonmin)
            comp,comp_point,count = 
	      composite_calc(time,comp_num,data,comp,comp_point,count)
          end
        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 composite_calc(time,comp_num,data,comp,comp_point,count)
    
    # 軸の値を保存
    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

    print '(lon,  time ) = (',comp_num, ',', time, ')', "\n"
    time = (time*4).to_i
    comp_num = comp_num*grid_lon_size/360 + grid_lon_size/2 
    print '(lgrid,tgrid) = (',comp_num, ',', time, ')', "\n"
    comp_point[comp_num,time] = 1.0
    
    # 作業配列
    tmp = NArray.sfloat(grid_lon_size,grid_sigmap_size).fill!(0.0)
    n =  ((grid_lon_size/2) - comp_num ).to_i
    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 

    return comp,comp_point,count

  end
    
end