#!/usr/bin/env ruby
# ----------------------------------------------
# local load path

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

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

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

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


END{


  setopt
  if ARGV.index("-marge")
    mknc_file_marge
  elsif ARGV.index("-horizontal-marge")
    mknc_file_marge_horizontal
  elsif ARGV.index("-horizontal")
    mknc_spctfilter_horizontal
  else
    mknc_spctfilter
  end

}

# ----------------------------------------------
# 引数解析
def setopt

  num = 0


  $cum = "adj"
  $expid = "con"
  $filter = "kelvin"  # "grav", "advect", "non"

  $cum = ARGV[num+1]      if num = ARGV.index("-cum")
  $expid = ARGV[num+1]      if num = ARGV.index("-exp")
  $var = ARGV[num+1]      if num = ARGV.index("-var")
  $filter = ARGV[num+1]      if num = ARGV.index("-fil")

  print "$cum = #{$cum}, $exp = #{$expid}, $var = #{$var}, #{Time.now}\n"

end

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


def mknc_file_marge

  $gr2ncfile_path = "./tmp/"
  file =  
    NetCDF.create("AGCM5-#{$cum}-#{$expid}_#{$filter}filt.nc")
  puts "AGCM5-#{$cum}-#{$expid}_#{$filter}filt.nc create start"

  file_name = Array.new
  Dir.foreach($gr2ncfile_path) { |item|
    file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /AGCM5-#{$cum}-#{$expid}_#{$filter}filt_/
  }

  file_name.each{ |moto|

    puts moto
    $t = Ape.new(moto)

    $t.netcdf_open.var_names.each { |name| 

      unless name == "x"  || name == "z"  || name == "time" || name == "lon" || name == "lat" || name == "sigmam"|| name == "sigmap" || name =~ /_wgt/ 
        puts name
        gphys = $t.go(name)
        
        # netCDF 初期化
        GPhys::NetCDF_IO.write(file, gphys )
      
      end
    }      

  }
    # 大域属性
    file.put_att("Conventions", "CF-1.0")
  file.put_att("title","Aqua Planet: AGCM5-#{$cum}-#{$expid} Experiment, #{$filter} filter")
    file.put_att("history", "Original data produced: #{Time.now} from TR netCDF by Yukiko YAMADA (GFD Dennou Club)")
    file.put_att("institution", "agcm5.3-20020620.tgz")
    file.close

end

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

def mknc_file_marge_horizontal

  $gr2ncfile_path = "./tmp/"
  file =  
    NetCDF.create("AGCM5-#{$cum}-#{$expid}_#{$filter}filt-horizontal.nc")
  puts "AGCM5-#{$cum}-#{$expid}_#{$filter}filt-horizontal.nc create start"

  file_name = Array.new
  Dir.foreach($gr2ncfile_path) { |item|
    file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /AGCM5-#{$cum}-#{$expid}_#{$filter}filt-horizontal_/
  }

  file_name.each{ |moto|
    
    puts moto
    $t = Ape.new(moto)

    $t.netcdf_open.var_names.each { |name| 

      unless name == "x"  || name == "z"  || name == "time" || name == "lon" || name == "lat" || name == "sigmam"|| name == "sigmap" || name =~ /_wgt/ 
        puts name
        gphys = $t.go(name)
        
        # netCDF 初期化
        GPhys::NetCDF_IO.write(file, gphys )
        
      end
    }      

  }
    # 大域属性
    file.put_att("Conventions", "CF-1.0")
  file.put_att("title","Aqua Planet: AGCM5-#{$cum}-#{$expid} Experiment, #{$filter} filter")
    file.put_att("history", "Original data produced: #{Time.now} from TR netCDF by Yukiko YAMADA (GFD Dennou Club)")
    file.put_att("institution", "agcm5.3-20020620.tgz")
    file.close
  
end

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

def mknc_spctfilter

  $var_array = [
    "tr_om", "tr_q", "tr_t", "tr_tconv", "tr_tppn", "tr_u", "tr_v"
  ]

  $var_gt3file_hash = { 
    "tr_om"    => "sigdot_out_y32" , 
    "tr_q"     => "q_out_y32", 
    "tr_t"     => "t_out_y32" , 
    "tr_tconv" => "qcnd_out_y32" , 
    "tr_tppn"  => "rain_sum_y32" , 
    "tr_u"     => "u_out_y32" , 
    "tr_v"     => "v_out_y32" 
  }

  $var_gt3filevar_hash = { 
    "tr_om"    => "sigdot", 
    "tr_q"     => "q", 
    "tr_t"     => "t", 
    "tr_tconv" => "qcnd", 
    "tr_tppn"  => "rain" , 
    "tr_u"     => "u", 
    "tr_v"     => "v"
  }
  
  $var_gt3fileunits_hash = { 
    "tr_om"    => "s-1", 
    "tr_q"     => "1", 
    "tr_t"     => "K", 
    "tr_tconv" => "K s-1", 
    "tr_tppn"  => "kg m-2 s-1" , 
    "tr_u"     => "m s-1", 
    "tr_v"     => "m s-1"
  }
  
  $var_gt3fileapename_hash = { 
    "tr_om"    => "sigma-vertical_vel.", 
    "tr_q"     => "specific_humidity", 
    "tr_t"     => "air_temperature", 
    "tr_tconv" => "condensation_heating", 
    "tr_tppn"  => "precipitation_flux" , 
    "tr_u"     => "eastward_wind", 
    "tr_v"     => "northward_wind"
  }

  file = NetCDF.create("AGCM5-#{$cum}-#{$expid}_#{$filter}filt_#{$var.downcase}.nc")
  puts "AGCM5-#{$cum}-#{$expid}_#{$filter}filt_#{$var.downcase}.nc create start"

  $gt3dir = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/#{$cum}700#{$expid}"

  infile = [
    "#{$gt3dir}400/#{$var_gt3file_hash[$var]}.nc",
    "#{$gt3dir}500/#{$var_gt3file_hash[$var]}.nc",
    "#{$gt3dir}600/#{$var_gt3file_hash[$var]}.nc",
    "#{$gt3dir}700/#{$var_gt3file_hash[$var]}.nc" 
    ]

  if $var == "tr_tppn"

    gphys = 
#      GPhys::NetCDF_IO.open(infile, $var_gt3filevar_hash[$var])[true,-1440..-1]
      GPhys::NetCDF_IO.open(infile, $var_gt3filevar_hash[$var])[true,0..1439]
    
    data = VArray.new(gphys.data.val.to_na*0.4*10**(-6)).rename($var)
    data = data.set_att("units", $var_gt3fileunits_hash[$var])
    data = data.set_att("ape_name",$var_gt3fileapename_hash[$var])
    
    g0 = VArray.new(gphys.grid_copy.axis(0).pos.val.to_na ).rename("lon")
    g0 = g0.set_att("standard_name","longitude")
    g0 = g0.set_att("long_name","longitude")
    g0 = g0.set_att("units","degrees_east")
    g0 = Axis.new().set_pos(g0)
    
    g2 = VArray.new( NArray.sfloat(1440).indgen!(1)*0.25 ).rename("time")
    g2 = g2.set_att("standard_name","time")
    g2 = g2.set_att("long_name","time")
    g2 = g2.set_att("units","days since 0000-01-01")
    g2 = g2.set_att("bounds","time_bnds")
    g2 = Axis.new().set_pos(g2)

    grid = Grid.new(g0,g2)

  else
    gphys = 
#      GPhys::NetCDF_IO.open(infile, $var_gt3filevar_hash[$var])[true,true,-1440..-1]
      GPhys::NetCDF_IO.open(infile, $var_gt3filevar_hash[$var])[true,true,0..1439]

    data = VArray.new(gphys.data.val.to_na).rename($var)
    data = data.set_att("units", $var_gt3fileunits_hash[$var])
    data = data.set_att("ape_name",$var_gt3fileapename_hash[$var])
    
    g0 = VArray.new(gphys.grid_copy.axis(0).pos.val.to_na ).rename("lon")
    g0 = g0.set_att("standard_name","longitude")
    g0 = g0.set_att("long_name","longitude")
    g0 = g0.set_att("units","degrees_east")
    g0 = Axis.new().set_pos(g0)
    
    g1 = VArray.new(gphys.grid_copy.axis(1).pos.val.to_na ).
      rename(gphys.grid_copy.axis(1).name)
    g1 = g1.set_att("standard_name","sigma")
    g1 = g1.set_att("long_name","sigma")
    g1 = g1.set_att("units","1")
    g1 = Axis.new().set_pos(g1)
        
    g2 = VArray.new( NArray.sfloat(1440).indgen!(0)*0.25 ).rename("time")
    g2 = g2.set_att("standard_name","time")
    g2 = g2.set_att("long_name","time")
    g2 = g2.set_att("units","days since 0000-01-01")
    g2 = g2.set_att("bounds","time_bnds")
    g2 = Axis.new().set_pos(g2)

    grid = Grid.new(g0,g1,g2)

  end
  
  gphys = GPhys.new(grid,data).set_lost_axes("lat=1.3953 degrees_north")

  if $filter =~ /^k/
    gphys_out = kelvin_filt(gphys)  
  elsif $filter =~ /^a/
    gphys_out = advect_filt(gphys)  
  elsif $filter =~ /^g/
    gphys_out = grav_filt(gphys)  
  elsif $filter =~ /^n/
    gphys_out = gphys[false,0..400]  
  end
    
  gphys_out = gphys_out.set_att("ape_name",gphys.get_att("ape_name")). 
                set_att("units",gphys.get_att("units"))

  if $var == "tr_tppn"
    
    data = VArray.new(gphys_out.data.val).rename(gphys_out.name)
    data = data.set_att("units", $var_gt3fileunits_hash[$var])
    data = data.set_att("ape_name",$var_gt3fileapename_hash[$var])
    
    g0 = VArray.new(gphys_out.grid_copy.axis(0).pos.val ).rename("lon")
    g0 = g0.set_att("standard_name","longitude")
    g0 = g0.set_att("long_name","longitude")
    g0 = g0.set_att("units","degrees_east")
    g0 = Axis.new().set_pos(g0)
    
    g2 = VArray.new( NArray.sfloat(401).indgen!(0)*0.25 ).rename("time")
    g2 = g2.set_att("standard_name","time")
    g2 = g2.set_att("long_name","time")
    g2 = g2.set_att("units","days since 0000-01-01")
    g2 = g2.set_att("bounds","time_bnds")
    g2 = Axis.new().set_pos(g2)

    grid = Grid.new(g0,g2)

  else
    
    data = VArray.new(gphys_out.data.val).rename(gphys_out.name)
    data = data.set_att("units", $var_gt3fileunits_hash[$var])
    data = data.set_att("ape_name",$var_gt3fileapename_hash[$var])
    
    g0 = VArray.new(gphys_out.grid_copy.axis(0).pos.val ).rename("lon")
    g0 = g0.set_att("standard_name","longitude")
    g0 = g0.set_att("long_name","longitude")
    g0 = g0.set_att("units","degrees_east")
    g0 = Axis.new().set_pos(g0)
    
    g1 = VArray.new(gphys_out.grid_copy.axis(1).pos.val ).
      rename(gphys_out.grid_copy.axis(1).name)
    g1 = g1.set_att("standard_name","sigma")
    g1 = g1.set_att("long_name","sigma")
    g1 = g1.set_att("units","1")
    g1 = Axis.new().set_pos(g1)
        
    g2 = VArray.new( NArray.sfloat(401).indgen!(1)*0.25 ).rename("time")
#    g2 = VArray.new(gphys_out.grid_copy.axis(2).pos.val ).rename("time")
    g2 = g2.set_att("standard_name","time")
    g2 = g2.set_att("long_name","time")
    g2 = g2.set_att("units","days since 0000-01-01")
    g2 = g2.set_att("bounds","time_bnds")
    g2 = Axis.new().set_pos(g2)

    grid = Grid.new(g0,g1,g2)

  end

  gphys_out = GPhys.new(grid,data).
    set_att("lost_axis","lat=1.3953 degrees_north")

#  if  gphys.lost_axes
#    gphys_out = gphys_out.
#      set_att("lost_axis", gphys.lost_axes.to_s)
#  end

  # netCDF 
  GPhys::NetCDF_IO.write(file, gphys_out )
#   GPhys::NetCDF_IO.write(file, spct )
  
  # 大域属性
  file.put_att("Conventions", "CF-1.0")
  file.put_att("title","Aqua Planet: #{$rezol}_#{$sstid} Experiment, #{$filter} filter")
  file.put_att("history", "Original data produced: #{Time.now} from TR netCDF by Yukiko YAMADA (AGU for APE)")
  file.put_att("institution", "#{$rezol}")
  file.close
  
end


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

def mknc_spctfilter_horizontal
  

  $var_array = [
    "tr_tppn", "tr_mslp", "tr_u_z12", "tr_u_z8", "tr_u_z5", "tr_u_z3", 
    "tr_v_z12", "tr_v_z8", "tr_v_z5", "tr_v_z3",
    "tr_gph_z12", "tr_gph_z5", "tr_t_z8", "tr_q_z5" 
  ]

  $var_gt3file_hash = { 
    "tr_tppn"    => "rain_sum",  
    "tr_mslp"    => "ps_out",  
    "tr_u_z12"   => "u_out_z12", 
    "tr_u_z8"    => "u_out_z8", 
    "tr_u_z5"    => "u_out_z5", 
    "tr_u_z3"    => "u_out_z3", 
    "tr_v_z12"   => "v_out_z12", 
    "tr_v_z8"    => "v_out_z8", 
    "tr_v_z5"    => "v_out_z5", 
    "tr_v_z3"    => "v_out_z3", 
    "tr_gph_z12" => "gph_out_z12", 
    "tr_gph_z5"  => "gph_out_z5", 
    "tr_t_z8"    => "t_out_z8", 
    "tr_q_z5"    => "q_out_z5"
  }

  $var_gt3filevar_hash = { 
    "tr_tppn"    => "rain",  
    "tr_mslp"    => "ps",  
    "tr_u_z12"   => "u", 
    "tr_u_z8"    => "u", 
    "tr_u_z5"    => "u", 
    "tr_u_z3"    => "u", 
    "tr_v_z12"   => "v", 
    "tr_v_z8"    => "v", 
    "tr_v_z5"    => "v", 
    "tr_v_z3"    => "v", 
    "tr_gph_z12" => "z", 
    "tr_gph_z5"  => "z", 
    "tr_t_z8"    => "t", 
    "tr_q_z5"    => "q"
  }
  
  $var_gt3fileunits_hash = { 
    "tr_tppn"    => "kg m-2 s-1",  
    "tr_mslp"    => "Pa",  
    "tr_u_z12"   => "m s-1", 
    "tr_u_z8"    => "m s-1", 
    "tr_u_z5"    => "m s-1", 
    "tr_u_z3"    => "m s-1", 
    "tr_v_z12"   => "m s-1", 
    "tr_v_z8"    => "m s-1", 
    "tr_v_z5"    => "m s-1", 
    "tr_v_z3"    => "m s-1", 
    "tr_gph_z12" => "m", 
    "tr_gph_z5"  => "m", 
    "tr_t_z8"    => "K", 
    "tr_q_z5"    => "1"
  }
  
  $var_gt3fileapename_hash = { 
    "tr_tppn"    => "precipitation_flux" , 
    "tr_mslp"    => "Pa",  
    "tr_u_z12"   => "eastward_wind", 
    "tr_u_z8"    => "eastward_wind", 
    "tr_u_z5"    => "eastward_wind", 
    "tr_u_z3"    => "eastward_wind", 
    "tr_v_z12"   => "northward_wind", 
    "tr_v_z8"    => "northward_wind", 
    "tr_v_z5"    => "northward_wind", 
    "tr_v_z3"    => "northward_wind", 
    "tr_gph_z12" => "geopotential_height",
    "tr_gph_z5"  => "geopotential_height",
    "tr_t_z8"    => "air_temperature",  
    "tr_q_z5"    => "specific_humidity"
  }

  $var_gt3filelostaxes_hash = { 
    "tr_tppn"    => nil , 
    "tr_mslp"    => nil ,  
    "tr_u_z12"   => "sigma=0.22953",
    "tr_u_z8"    => "sigma=0.54945", 
    "tr_u_z5"    => "sigma=0.82977", 
    "tr_u_z3"    => "sigma=0.94994", 
    "tr_v_z12"   => "sigma=0.22953",
    "tr_v_z8"    => "sigma=0.54945",  
    "tr_v_z5"    => "sigma=0.82977",  
    "tr_v_z3"    => "sigma=0.94994",  
    "tr_gph_z12" => "sigma=0.22953",
    "tr_gph_z5"  => "sigma=0.82977", 
    "tr_t_z8"    => "sigma=0.54945",    
    "tr_q_z5"    => "sigma=0.82977"
  }


  file = NetCDF.create("AGCM5-#{$cum}-#{$expid}_#{$filter}filt-horizontal_#{$var.downcase}.nc")
  puts "AGCM5-#{$cum}-#{$expid}_#{$filter}filt-horizontal_#{$var.downcase}.nc create start"

  $gt3dir = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/#{$cum}700#{$expid}"



  infile = [
    "#{$gt3dir}400/#{$var_gt3file_hash[$var]}.nc",
    "#{$gt3dir}500/#{$var_gt3file_hash[$var]}.nc",
    "#{$gt3dir}600/#{$var_gt3file_hash[$var]}.nc",
    "#{$gt3dir}700/#{$var_gt3file_hash[$var]}.nc" 
    ]
  gphys = 
    GPhys::NetCDF_IO.open(infile, $var_gt3filevar_hash[$var])[true,true,0..1439]

  if $var == "tr_tppn"
    data = VArray.new(gphys.data.val.to_na*0.4*10**(-6)).rename($var)
  elsif  $var == "tr_mslp"
    data = VArray.new(gphys.data.val.to_na*100.0).rename($var)
  else
    data = VArray.new(gphys.data.val.to_na).rename($var)
  end

  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units", $var_gt3fileunits_hash[$var])
  data = data.set_att("ape_name",$var_gt3fileapename_hash[$var])
  
  g0 = VArray.new(gphys.grid_copy.axis(0).pos.val.to_na ).rename("lon")
  g0 = g0.set_att("standard_name","longitude")
  g0 = g0.set_att("long_name","longitude")
  g0 = g0.set_att("units","degrees_east")
  g0 = Axis.new().set_pos(g0)
  
  g1 = VArray.new(gphys.grid_copy.axis(1).pos.val.to_na ).rename("lat")
  g1 = g1.set_att("standard_name","latitude")
  g1 = g1.set_att("long_name","latitude")
  g1 = g1.set_att("units","degrees_north")
  g1[-1..0] = g1[0..-1] # 緯度軸の反転
  g1 = Axis.new().set_pos(g1)

  g2 = VArray.new( NArray.sfloat(1440).indgen!(0)*0.25 ).rename("time")
  g2 = g2.set_att("standard_name","time")
  g2 = g2.set_att("long_name","time")
  g2 = g2.set_att("units","days since 0000-01-01")
  g2 = g2.set_att("bounds","time_bnds")
  g2 = Axis.new().set_pos(g2)

  grid = Grid.new(g0,g1,g2)
  gphys = GPhys.new(grid,data).set_lost_axes($var_gt3filelostaxes_hash[$var]).
    cut(true,-30..30,true)
  
  if $filter =~ /^k/
    gphys_out = kelvin_filt(gphys)  
  elsif $filter =~ /^a/
    gphys_out = advect_filt(gphys)  
  elsif $filter =~ /^g/
    gphys_out = grav_filt(gphys)  
  elsif $filter =~ /^n/
    gphys_out = gphys[false,0..400]  

  end
    
  gphys_out = gphys_out.set_att("ape_name",gphys.get_att("ape_name")). 
                set_att("units",gphys.get_att("units"))

  data = VArray.new(gphys_out.data.val).rename(gphys_out.name)
  data = data.set_att("units", $var_gt3fileunits_hash[$var])
  data = data.set_att("ape_name",$var_gt3fileapename_hash[$var])
  
  g0 = VArray.new(gphys_out.grid_copy.axis(0).pos.val ).rename("lon")
  g0 = g0.set_att("standard_name","longitude")
  g0 = g0.set_att("long_name","longitude")
  g0 = g0.set_att("units","degrees_east")
  g0 = Axis.new().set_pos(g0)
  
  g1 = VArray.new(gphys.grid_copy.axis(1).pos.val ).rename("lat")
  g1 = g1.set_att("standard_name","latitude")
  g1 = g1.set_att("long_name","latitude")
  g1 = g1.set_att("units","degrees_north")
  g1 = Axis.new().set_pos(g1)
        
  g2 = VArray.new( NArray.sfloat(401).indgen!(0)*0.25 ).rename("time")
#  g2 = VArray.new(gphys_out.grid_copy.axis(2).pos.val ).rename("time")
  g2 = g2.set_att("standard_name","time")
  g2 = g2.set_att("long_name","time")
  g2 = g2.set_att("units","days since 0000-01-01")
  g2 = g2.set_att("bounds","time_bnds")
  g2 = Axis.new().set_pos(g2)
  
  grid = Grid.new(g0,g1,g2)

  gphys_out = GPhys.new(grid,data)

  if $var == "tr_tppn" || $var == "tr_mslp"
  else
    gphys_out = gphys_out.
      set_att("lost_axis", $var_gt3filelostaxes_hash[$var])
  end

  # netCDF 
  GPhys::NetCDF_IO.write(file, gphys_out )
  # GPhys::NetCDF_IO.write(file, spct )
  
  # 大域属性
  file.put_att("Conventions", "CF-1.0")
  file.put_att("title","Aqua Planet: #{$rezol}_#{$sstid} Experiment, #{$filter} filter")
  file.put_att("history", "Original data produced: #{Time.now} from TR netCDF by Yukiko YAMADA (AGU for APE)")
  file.put_att("institution", "#{$rezol}")
  file.close
  
end

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

def advect_filt(gphys)  

#  spct = gphys.fft(nil,0,-1)
#  spct_tmp = spct.copy

# $x_size = spct.grid_copy.axis(0).pos.val  

  if gphys.rank == 2

    spct = gphys.fft(nil,0,-1)
    spct_tmp = spct.copy
    $x_size = spct.grid_copy.axis(0).pos.val  
    spct = spct*0.0

    (4..28).each{ |num|
      
      up = 360/25 * 5.1/30 * num
      bo = 360/25 * 1.01/30 * num

      up = 0.51*360.0/4 if up >  0.51*360.0/4

      work = 
        (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
        spct_tmp.cut(num,false,bo..up).data.to_a +
        (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
      work = NArray.to_na(work)
      
      spct[num,false,true] = spct[num,false,true] + work    
    }

    ( ($x_size.size-28)..($x_size.size-4) ).each{ |num|
      
      up =  ( 360/25 * 1.01/30 * ( - $x_size.size + num ) ) +360
      bo =  ( 360/25 * 5.1/30 * ( - $x_size.size + num ) ) +360

      bo = 360 - 0.51*360.0/4 if bo < ( 360 - 0.51*360.0/4 )

      work = 
        (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
        spct_tmp.cut(num,false,bo..up).data.to_a +
        (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
      work = NArray.to_na(work)
      
      spct[num,false,true] = spct[num,false,true] + work    
    }

    gphys = spct.fft(true).real

    gphys = gphys[true,0..400]
  else

    grid_0 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(gphys.grid_copy.axis(0).pos.val.size).indgen!).rename("lon").
      put_att("units","knum"))
    grid_1 = gphys.grid_copy.axis(1)	
    grid_2 = gphys.grid_copy.axis(2)	
    grid = Grid.new(grid_0,grid_1,grid_2)
    data = VArray.new(NArray.complex(grid_0.pos.val.size,grid_1.pos.val.size,grid_2.pos.val.size).fill!(0.0)).rename(gphys.name)
    spct = GPhys.new(grid, data )

    levsize = spct.grid_copy.axis(1).pos.val.size
    levsize.times{ |levnum|

      print "lev=#{levnum}, #{Time.now}\n"
      
       spct[true,levnum,true] = gphys[true,levnum,true].fft(nil)
       spct_tmp = spct[true,levnum,true].copy
       $x_size = spct.grid_copy.axis(0).pos.val  
       spct[true,levnum,true] = 0.0

      (4..28).each{ |num|
        
        up = 360/25 * 5.1/30 * num
        bo = 360/25 * 1.01/30 * num
        
        up = 0.51*360.0/4 if up >  0.51*360.0/4

        work = 
          (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up).data.to_a +
          (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
        work = NArray.to_na(work)
        
        spct[num,levnum,true] = spct[num,levnum,true] + work    
      }

      ( ($x_size.size-28)..($x_size.size-4) ).each{ |num|
        
        up = ( 360.0/25 * 1.01/30 * (- $x_size.size*1.0 + num ) ) +360
        bo = ( 360.0/25 * 5.1/30 * (- $x_size.size*1.0 + num ) ) +360
        
      bo = 360 - 0.51*360.0/4 if bo < ( 360 - 0.51*360.0/4 )

        work = 
          (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up).data.to_a +
          (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
        work = NArray.to_na(work)
        
        spct[num,levnum,true] = spct[num,levnum,true] + work    
      }

      gphys[true,levnum,true] = spct[true,levnum,true].fft(true).real
      
    }

      gphys = gphys[true,true,0..400]

  end

  gphys = gphys.rename("#{gphys.name}_adfilt")
  return gphys
  
end

def kelvin_filt(gphys)  


  if gphys.rank == 2

    spct = gphys.fft(nil,0,-1)
    spct_tmp = spct.copy
    $x_size = spct.grid_copy.axis(0).pos.val  
    spct = spct*0.0

    kelvin_bo, grav1 = bunsan_calc(8)
  #  kelvin_up, grav1 = bunsan_calc(90)
    kelvin_up, grav1 = bunsan_calc(200)

    ( ($x_size.size-15)..($x_size.size-3) ).each{ |num|
      
      up = kelvin_up.data.val[$x_size.size-num] * 360.0/4
      bo = kelvin_bo.data.val[$x_size.size-num] * 360.0/4
      
      up = 0.41*360.0/4 if up >  0.41*360.0/4
      
      work = 
        (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
        spct_tmp.cut(num,false,bo..up).data.to_a +
        (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
      work = NArray.to_na(work)
      
      spct[num,false,true] = spct[num,false,true] + work     
    }
    
    ( 3..15 ).each{ |num|

      up = - kelvin_bo.data.val[num] * 360.0/4 + 360
      bo = - kelvin_up.data.val[num] * 360.0/4 + 360
      
      bo = 360- 0.41*360.0/4 if bo <  (360-0.41*360.0/4 )

      work = 
        (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
        spct_tmp.cut(num,false,bo..up).data.to_a +
        (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
      work = NArray.to_na(work)
      
      spct[num,false,true] = spct[num,false,true] + work     
    }


    gphys = spct.fft(true).real
    gphys = gphys[true,0..400]

  else

    grid_0 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(gphys.grid_copy.axis(0).pos.val.size).indgen!).rename("lon").
      put_att("units","knum"))
    grid_1 = gphys.grid_copy.axis(1)	
    grid_2 = gphys.grid_copy.axis(2)	
    grid = Grid.new(grid_0,grid_1,grid_2)
    data = VArray.new(NArray.complex(grid_0.pos.val.size,grid_1.pos.val.size,grid_2.pos.val.size).fill!(0.0)).rename(gphys.name)
    spct = GPhys.new(grid, data )

    $x_size = spct.grid_copy.axis(0).pos.val  
    kelvin_bo, grav1 = bunsan_calc(8)
  #  kelvin_up, grav1 = bunsan_calc(90)
    kelvin_up, grav1 = bunsan_calc(200)

    levsize = spct.grid_copy.axis(1).pos.val.size
    levsize.times{ |levnum|

      print "lev=#{levnum}, #{Time.now}\n"
      
       spct[true,levnum,true] = gphys[true,levnum,true].fft(nil)
       spct_tmp = spct[true,levnum,true].copy
       $x_size = spct.grid_copy.axis(0).pos.val  
       spct[true,levnum,true] = 0.0

      ( ($x_size.size-15)..($x_size.size-3) ).each{ |num|
  
        up = kelvin_up.data.val[$x_size.size-num] * 360.0/4
        bo = kelvin_bo.data.val[$x_size.size-num] * 360.0/4
      
        up = 0.41*360.0/4 if up >  0.41*360.0/4
        
        work = 
          (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up).data.to_a +
          (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
        work = NArray.to_na(work)
        
        spct[num,levnum,true] = spct[num,levnum,true] + work    
  
       }

      ( 3..15 ).each{ |num|
  
        up = - kelvin_bo.data.val[num] * 360.0/4 + 360
        bo = - kelvin_up.data.val[num] * 360.0/4 + 360
            
        bo = 360- 0.41*360.0/4 if bo <  (360 - 0.41*360.0/4 )
       
        work = 
          (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up).data.to_a +
          (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
        work = NArray.to_na(work)
        
        spct[num,levnum,true] = spct[num,levnum,true] + work    
  
       }

      gphys[true,levnum,true] = spct[true,levnum,true].fft(true).real
      
    }

    gphys = gphys[true,true,0..400]

  end
  
  gphys = gphys.rename("#{gphys.name}_kfilt")
  return gphys

end


def grav_filt(gphys)  

  if gphys.rank == 2

    spct = gphys.fft(nil,0,-1)
    spct_tmp = spct.copy
    $x_size = spct.grid_copy.axis(0).pos.val  
    spct = spct*0.0

    $x_size = spct.grid_copy.axis(0).pos.val * (-1.0)
    kelvin, grav1_bo = bunsan_calc(12)
    kelvin, grav1_up = bunsan_calc(50)

    (1..15).each{ |num|
      
      up = grav1_up.data.val[num] * 360.0/4
      bo = grav1_bo.data.val[num] * 360.0/4
      
      up = 0.71*360.0/4 if up >  0.71*360.0/4 

      work = 
        (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
        spct_tmp.cut(num,false,bo..up).data.to_a +
        (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
      work = NArray.to_na(work)
      
      spct[num,false,true] = spct[num,false,true] + work    
      
    }

    ( ($x_size.size-15)..($x_size.size-1) ).each{ |num|

      up = - grav1_bo.data.val[$x_size.size-num] * 360.0/4 +360
      bo = - grav1_up.data.val[$x_size.size-num] * 360.0/4 +360
      
      bo = 360 -0.71*360.0/4 if bo <  (360 -0.71*360.0/4 )

      work = 
        (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
        spct_tmp.cut(num,false,bo..up).data.to_a +
        (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
      work = NArray.to_na(work)
      
      spct[num,false,true] = spct[num,false,true] + work    
      
    }

    gphys = spct.fft(true).real
    gphys = gphys[true,0..400]

  else

      grid_0 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(gphys.grid_copy.axis(0).pos.val.size).indgen!).rename("lon").
      put_att("units","knum"))
    grid_1 = gphys.grid_copy.axis(1)	
    grid_2 = gphys.grid_copy.axis(2)	
    grid = Grid.new(grid_0,grid_1,grid_2)
    data = VArray.new(NArray.complex(grid_0.pos.val.size,grid_1.pos.val.size,grid_2.pos.val.size).fill!(0.0)).rename(gphys.name)
    spct = GPhys.new(grid, data )

    $x_size = spct.grid_copy.axis(0).pos.val * (-1.0)
    kelvin, grav1_bo = bunsan_calc(12)
    kelvin, grav1_up = bunsan_calc(50)

    levsize = spct.grid_copy.axis(1).pos.val.size
    levsize.times{ |levnum|

      print "lev=#{levnum}, #{Time.now}\n"
      
       spct[true,levnum,true] = gphys[true,levnum,true].fft(nil)
       spct_tmp = spct[true,levnum,true].copy
       $x_size = spct.grid_copy.axis(0).pos.val  
       spct[true,levnum,true] = 0.0

      (1..15).each{ |num|
    
        up = grav1_up.data.val[num] * 360.0/4
        bo = grav1_bo.data.val[num] * 360.0/4

        up = 0.71*360.0/4 if up >  0.71*360.0/4 

        work = 
          (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up).data.to_a +
          (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
        work = NArray.to_na(work)

        spct[num,levnum,true] = spct[num,levnum,true] + work    
      }

      ( ($x_size.size-15)..($x_size.size-1) ).each{ |num|

        up = - grav1_bo.data.val[$x_size.size-num] * 360.0/4 +360
        bo = - grav1_up.data.val[$x_size.size-num] * 360.0/4 +360
      
        bo = 360 -0.71*360.0/4 if bo <  (360 -0.71*360.0/4 )

        work = 
          (spct_tmp.cut(num,false,0..bo).data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up).data.to_a +
          (spct_tmp.cut(num,false,up..999).data*0.0).to_a 
        work = NArray.to_na(work)

        spct[num,levnum,true] = spct[num,levnum,true] + work    
      }

      gphys[true,levnum,true] = spct[true,levnum,true].fft(true).real
      
    }

    gphys = gphys[true,true,0..400]

  end
  
  gphys = gphys.rename("#{gphys.name}_gfilt")
  return gphys

end

def bunsan_calc(h=200)
  
  omega = 2*3.1415/24/60/60
  beta =  2*(omega)/(4e+7) * ( 60*60*24*(6370e+3) )  # (2*Omega/a)*cos(0)
  c_g = ((h*9.8)**0.5)*60*60*24/(4e+7)
  
  # 軸
  knum = VArray.new($x_size).rename("knum").put_att("units","1")
  knum = Axis.new().set_pos(knum)
  grid = Grid.new(knum)
        
  ## 虚数を含む軸に変換
  x_size = NArray.complex($x_size.size).fill!(1.0) * $x_size

  # モード
  $sym_num = 1
  num = $sym_num
    
  ## 3 次方程式解法: カルダノの公式
  ## x^3 + 3px + q
  ## t^2 + qt - p =0, t = -q/2 +- (q**2/4 + p)**0.5
  p = - ( (c_g*x_size)**2 + c_g*beta*(2*num+1)) / 3
  q = - c_g**2*beta*x_size
  t1  =  -q/2 + (q**2/4 + p**3)**0.5
  t2  =  -q/2 - (q**2/4 + p**3)**0.5
  t1 = t1**(1.0/3)
  t2 = t2**(1.0/3)
  
  # 重力波
  $grav1 = (t1 + t2).real
    $grav1 = VArray.new($grav1).rename("grav1").
    put_att("units","1").put_att("ape_name","gravity")
  $grav1 = GPhys.new( grid, $grav1)
  
  # 反対重力波 (図示しない)
  $grav2 = (t1 *  Complex.polar(1, Math::PI*2/3) + 
	     t2 * Complex.polar(1, Math::PI*4/3)).real
  $grav2 = VArray.new($grav2).rename("grav2").
    put_att("units","1").put_att("ape_name","gravity")
  $grav2 = GPhys.new( grid, $grav2)
  
  # ロスビー波
  $rossby = (t2 * Complex.polar(1, Math::PI*2/3) + 
	      t1 * Complex.polar(1, Math::PI*4/3) ).real
  $rossby2 = $rossby
  $rossby = ($rossby > 0) * $rossby
    $rossby = VArray.new($rossby).rename("rossby").
    put_att("units","1").put_att("ape_name","rossby")
  $rossby = GPhys.new( grid, $rossby)
  
  $rossby2 = VArray.new($rossby2).rename("rossby").
    put_att("units","1").put_att("ape_name","rossby")
  $rossby2 = GPhys.new( grid, $rossby2)
  
  # ケルビン波
  $kelvin = (c_g*x_size).real
  $kelvin = ($kelvin > 0) * $kelvin
  $kelvin = VArray.new($kelvin).rename("kelvin").
    put_att("units","1").put_att("ape_name","kelvin")
  $kelvin = GPhys.new( grid, $kelvin)
  
  # 混合ロスビー
  $mixed = c_g*x_size/2 + ((c_g*x_size/2)**2 + c_g*beta )**0.5
  $mixed = VArray.new($mixed).rename("mixed").
    put_att("units","1").put_att("ape_name","mixed")
  $mixed = GPhys.new( grid, $mixed)

  return $kelvin , $grav1

end