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

 $local_path = '/home/yukiko/work/ape/yukiko/lib'
# $local_path = '/home/yukiko/tmp/ape-data/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

  $sstid = "control" 
  $rezol = "T159L48_eml"
  $filter = "kelvin"  # "grav", "advect"

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

  $groupid = $groupid_hash[$rezol]

  $num = $num.to_i
end

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


def mknc_file_marge

  $gr2ncfile_path = "./tmp/"
  file =  
    NetCDF.create("#{$rezol}_#{$sstid}_#{$filter}filt.nc")
  puts "#{$rezol}_#{$sstid}_#{$filter}filt.nc create start"

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

  file_name.each{ |moto|

    puts moto
    $t = Ape.new(moto)

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

      if name =~ /tr_/
          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: #{$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_file_marge_horizontal

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

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

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

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

      if name =~ /tr_/
        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: #{$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

  name_array = [
    [00,"tr_tppn"  , "precipitation_flux","kg m-2 s-1"], 
    [11,"tr_u"     , "eastward_wind","m s-1"], 
    [12,"tr_v"     , "northward_wind","m s-1"], 
    [13,"tr_t"     , "air_temperature","K"], 
    [14,"tr_om"    , "omega","Pa s-1"], 
    [16,"tr_q"     , "specific_humidity","1"], 
    [22,"tr_t_conv", "Temperature tendency - Convection","K s-1"], 
    [23,"tr_t_cld" , "Temperature tendency - Cloud","K s-1"] 
    ]
  
  $num = $num.to_i

  print "$rezol = #{$rezol}, #{name_array[$num][1]}, #{Time.now}\n"

  file = NetCDF.create("#{$rezol}_#{$sstid}_#{$filter}filt_#{name_array[$num][1]}.nc")
  puts "#{$rezol}_#{$sstid}_#{$filter}filt_#{name_array[$num][1]}.nc create start"

  $gr2ncfile_path = "/home/yukiko/tr_request/#{$rezol}/"

  if $num == 0 
    @data = Ape.new("#{$gr2ncfile_path}/#{$groupid}_TR00_#{$sstid}.nc") 
  else
    @data = Ape.new("#{$gr2ncfile_path}/#{$groupid}_TR#{name_array[$num][0]}_#{$sstid}.nc") 
  end
#    @data = Ape.new("#{$rezol}_#{$sstid}_nofilt_#{name_array[$num][1]}.nc")
    
  if $rezol == "ECMWF"
    time = 1459
  elsif $rezol == "GSFC" || $rezol == "LASG"
    time = 1460
  elsif $rezol == "NCAR"
    time = 1461
  elsif $rezol == "AGUforAPE"
    time = 1440
  end
  
  if name_array[$num][1] == "tr_tppn"
    gphys = @data.gphys_open(name_array[$num][1]).cut(true,0,true)
    grid_1 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(time).indgen!/4).rename("time").
              put_att("units","days")) 
    grid_0 = gphys.grid_copy.axis(0)	
    grid = Grid.new(grid_0,grid_1)
    gphys = GPhys.new(grid, gphys.data ).set_lost_axes(gphys.lost_axes)
  else

    gphys = @data.gphys_open(name_array[$num][1]).cut(true,0,true,true)

##   if name_array[$num][1] == "tr_t_cld" || name_array[$num][1] == "tr_t_conv"
#     ps = Ape.new("#{$gr2ncfile_path}/#{$groupid}_TR00_#{$sstid}.nc").
#                         gphys_open("tr_mslp").cut(true,0,true)
#     plev = NArray[1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000, 25000, 30000, 40000, 50000, 60000, 70000, 85000, 92500, 100000]
#     gphys = mks2p(gphys,ps,plev)
##   end

    grid_1 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(time).indgen!/4).rename("time").
              put_att("units","days")) 
    grid_0 = gphys.grid_copy.axis(0)	
    grid_00 = gphys.grid_copy.axis(1)	
    grid = Grid.new(grid_0,grid_00,grid_1)
    gphys = GPhys.new(grid, gphys.data ).set_lost_axes(gphys.lost_axes)
  end
#
#  gphys = @data.gphys_open("#{name_array[$num][1]}_nfilt")

  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/
#    if name_array[$num][1] == "tr_tppn"
#      gphys = gphys[true,-401..-1]
#    else
#      gphys = gphys[true,true,-401..-1]
#    end
    gphys_out = gphys.rename("#{gphys.name}_nfilt")
  end

  gphys_out = gphys_out.rename(gphys_out.name.gsub("_nfilt",""))

  #  gphys = spct.fft(true,0,-1).real.rename(spct.name.gsub("_spct","") ) 
  #  spct = spct.abs ** 2
    
#

#  gphy_out = gpyhs_out[true,true,-401..-1]

  gphys_out = gphys_out.set_att("ape_name",name_array[$num][2] ). 
                set_att("units",name_array[$num][3] )

  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
    
  name_array = [
#    [00,"tr_tppn", 0,"tr_tppn", "precipitation_flux","kg m-2 s-1"], 
#    [00,"tr_mslp", 0, "tr_mslp", "air_pressure_at_sea_level","Pa"], 
    [01,"tr_tppn", 0,"tr_tppn", "precipitation_flux","kg m-2 s-1"], 
    [06,"tr_mslp", 0, "tr_mslp", "air_pressure_at_sea_level","Pa"], 
    [13,"tr_t",60000,"tr_t600","air_temperature","K"], 
    [16,"tr_q",85000,"tr_q850","specific_humidity","1"], 
#--
    [11,"tr_u",25000,"tr_u250","eastward_wind","m s-1"], 
    [11,"tr_u",85000,"tr_u850","eastward_wind","m s-1"], 
    [11,"tr_u",60000,"tr_u600","eastward_wind","m s-1"], 
    [11,"tr_u",92500,"tr_u925","eastward_wind","m s-1"], 
#---
    [12,"tr_v",25000,"tr_v250","northward_wind","m s-1"], 
    [12,"tr_v",85000,"tr_v850","northward_wind","m s-1"], 
    [12,"tr_v",60000,"tr_v600","northward_wind","m s-1"], 
    [12,"tr_v",92500,"tr_v925","northward_wind","m s-1"], 
#---
    [15,"tr_z",25000,"tr_phi250","geopotential_height","m"], 
    [15,"tr_z",85000,"tr_phi850","geopotential_height","m"]
  ]
    
  $num = $num.to_i

  file = NetCDF.create("#{$rezol}_#{$sstid}_#{$filter}filt-horizontal_#{name_array[$num][3]}.nc")
  puts "#{$rezol}_#{$sstid}_#{$filter}filt-horizontal_#{name_array[$num][3]}.nc create start"

  $gr2ncfile_path = "/home/yukiko/tr_request/#{$rezol}/"

  if name_array[$num][1] == "tr_tppn" || name_array[$num][1] == "tr_mslp"
    if $rezol == "CSIRO"
      @data = Ape.new("#{$gr2ncfile_path}/#{$groupid}_TR0#{name_array[$num][0]}_#{$sstid}.nc") 
    else
      @data = Ape.new("#{$gr2ncfile_path}/#{$groupid}_TR00_#{$sstid}.nc") 
    end
  else
    @data = Ape.new("#{$gr2ncfile_path}/#{$groupid}_TR#{name_array[$num][0]}_#{$sstid}.nc") 
  end

  if name_array[$num][1] == "tr_tppn" || name_array[$num][1] == "tr_mslp"
    gphys = @data.gphys_open(name_array[$num][1])
  else
    lev = name_array[$num][2]
    lev = lev /100    if $rezol =~ /ECMWF/ || $rezol == "LASG"|| $rezol == "CSIRO"
    opn_name = name_array[$num][1]
    opn_name = "tr_phi" if $rezol == "AGUforAPE" && opn_name == "tr_z"
    gphys = @data.gphys_open(opn_name).cut(true,true,lev,true)
  end

    
  if $rezol =~ /ECMWF/
    time = 1459
  elsif $rezol == "GSFC" || $rezol == "LASG"|| $rezol == "CSIRO"
    time = 1460
  elsif $rezol == "NCAR"
    time = 1461
  elsif $rezol == "AGUforAPE"
    time = 1440
  end

  grid_1 =
    Axis.new().
    set_pos(VArray.new(NArray.sfloat(time).indgen!/4).rename("time").
            put_att("units","days")) 
  grid_0 = gphys.grid_copy.axis(0)	
  grid_00 = gphys.grid_copy.axis(1)	
  grid = Grid.new(grid_0,grid_00,grid_1)
  gphys_tmp = VArray.new( gphys.data.val.to_na).rename(gphys.name) 
  gphys = GPhys.new(grid, gphys_tmp ).set_lost_axes(gphys.lost_axes)

  gphys_out = gphys[true,true,-401..-1]*0.0

  if $filter =~ /^k/
    grid_00.pos.val.size.times{ |num|
    p "lev=#{num}"
    gphys_out[true,num,true] = kelvin_filt(gphys[true,num,true])  
  } 
   gphys_out = gphys_out.rename("#{name_array[$num][3]}_kfilt")
  elsif $filter =~ /^a/
    grid_00.pos.val.size.times{ |num|
    p "lev=#{num}"
    gphys_out[true,num,true] = advect_filt(gphys[true,num,true])  
  } 
    gphys_out = gphys_out.rename("#{name_array[$num][3]}_adfilt")
  elsif $filter =~ /^g/
    grid_00.pos.val.size.times{ |num|
    p "lev=#{num}"
    gphys_out[true,num,true] = grav_filt(gphys[true,num,true])  
  } 
    gphys_out = gphys_out.rename("#{name_array[$num][3]}_gfilt")
  elsif $filter =~ /^n/
    gphys_out = gphys[true,true,-401..-1]
    gphys_out = gphys_out.rename("#{name_array[$num][3]}_nfilt")
  end

  gphys_out = gphys_out.
  set_att("ape_name",name_array[$num][4] ).
  set_att("units",name_array[$num][5] )

  if gphys.lost_axes
    unless gphys.lost_axes.to_s == ""
    gphys_out = gphys_out.
      set_att("lost_axis", gphys.lost_axes.to_s)
    end
  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  

  gphys_tmp = gphys.copy

  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..15).each{ |num|
    (4..28).each{ |num|
      
      up = 360/25 * 5.1/30 * num
#      bo = 360/25 * 1.35/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-15)..($x_size.size-4) ).each{ |num|
    ( ($x_size.size-28)..($x_size.size-4) ).each{ |num|
      
#      up =  ( 360/25 * 1.35/30 * ( - $x_size.size + num ) ) +360
      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_tmp = spct.fft(true).real
  gphys_tmp = gphys_tmp[true,-401..-1]

  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"
      

#     if $rezol == "ECMWF" 
#       spct = gphys.to_na.fft(nil,0,-1)
#     else
#       spct = gphys.fft(nil,0,-1)
#     end
       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..15).each{ |num|
      (4..28).each{ |num|
        
        up = 360/25 * 5.1/30 * num
#        bo = 360/25 * 1.35/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-15)..($x_size.size-4) ).each{ |num|
      ( ($x_size.size-28)..($x_size.size-4) ).each{ |num|
        
#        up = ( 360.0/25 * 1.35/30 * (- $x_size.size*1.0 + num ) ) +360
        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_tmp[true,levnum,true] = spct[true,levnum,true].fft(true).real
  gphys_tmp = gphys_tmp[true,true,-401..-1]
      
    }

  end

#  spct = spct.rename("#{gphys.name}_spct_adfilt") 
#  gphys = spct.fft(true,0,-1).real.rename(spct.name.gsub("_spct","") ) 


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

def kelvin_filt(gphys)  

  gphys_tmp = gphys.copy

  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_bo, grav1 = bunsan_calc(2)
  #  kelvin_up, grav1 = bunsan_calc(90)
    kelvin_up, grav1 = bunsan_calc(200)

#    ( ($x_size.size-15)..($x_size.size-3) ).each{ |num|
    ( ($x_size.size-15)..($x_size.size-1) ).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|
    ( 1..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_tmp = spct.fft(true).real
  gphys_tmp = gphys_tmp[true,-401..-1]
  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_bo, grav1 = bunsan_calc(2)
  #  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|
      ( ($x_size.size-15)..($x_size.size-1) ).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|
      ( 1..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_tmp[true,levnum,true] = spct[true,levnum,true].fft(true).real
  gphys_tmp = gphys_tmp[true,true,-401..-1]  
      
    }

  end


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

end


def grav_filt(gphys)  

  gphys_tmp = gphys.copy

  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_tmp = spct.fft(true).real
  gphys_tmp = gphys_tmp[true,-401..-1]  
  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_tmp[true,levnum,true] = spct[true,levnum,true].fft(true).real
  gphys_tmp = gphys_tmp[true,true,-401..-1]        
    }

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

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_mslp
  g = @data.gphys_open("tr_mslp").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").sub("y=","lat="),
    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 __gr_tr_horizontal(var)
  g = @data.gphys_open(var)
  if g.data.get_att("lost_axis")
    lost_axis = [g.data.get_att("lost_axis").sub("y=","lat=")]
    g = g.set_lost_axes(lost_axis)
  end
  gphys = __gr2_gphys_make_3d_horizontal(g)
  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 $rezol =~ /T159L48/

    data[true,0..(3*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0011.nc",
                          g.name).cut(true,dimcut,true).data.val.to_na
    data[true,(3*30*4)..(3*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0012.nc",
                          g.name).cut(true,dimcut,true).data.val.to_na
    data[true,(3*30*4*2)..(3*30*4*3-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0013.nc",
                          g.name).cut(true,dimcut,true).data.val.to_na
    data[true,(3*30*4*3)..(3*30*4*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0014.nc",
                          g.name).cut(true,dimcut,true).data.val.to_na
  else

    data[true,0..(6*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0006.nc",
                          g.name).cut(true,dimcut,true).data.val.to_na
    data[true,(6*30*4)..(6*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0007.nc",
                          g.name).cut(true,dimcut,true).data.val.to_na
  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)

  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_0011.nc",
                          g.name).data.val.to_na
    data[true,true,(3*30*4)..(3*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0012.nc",
                          g.name).data.val.to_na
    data[true,true,(3*30*4*2)..(3*30*4*3-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0013.nc",
                          g.name).data.val.to_na
    data[true,true,(3*30*4*3)..(3*30*4*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0014.nc",
                          g.name).data.val.to_na
  else

    data[true,true,0..(6*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0006.nc",
                          g.name).data.val.to_na
    data[true,true,(6*30*4)..(6*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_0007.nc",
                          g.name).data.val.to_na
  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_horizontal(g)

  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_horizontal_0011.nc",
                          g.name).data.val.to_na
    data[true,true,(3*30*4)..(3*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_horizontal_0012.nc",
                          g.name).data.val.to_na
    data[true,true,(3*30*4*2)..(3*30*4*3-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_horizontal_0013.nc",
                          g.name).data.val.to_na
    data[true,true,(3*30*4*3)..(3*30*4*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_horizontal_0014.nc",
                          g.name).data.val.to_na
  else

    data[true,true,0..(6*30*4-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_horizontal_0006.nc",
                          g.name).data.val.to_na
    data[true,true,(6*30*4)..(6*30*4*2-1)] =
      GPhys::NetCDF_IO.open(
                          "#{$gr2ncfile_path}#{$rezol}_expID01_horizontal_0007.nc",
                          g.name).data.val.to_na
  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 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



def mks2p(gphys,ps,plev)

#  garr1 , gphys(im,jm,km)
#  garr2 , gphys_out(im,jm,plev.size)
#  gp1 , plev_org(im,jm,km)
#  gp2 , plev(plev.size)

  plev_org = sigma_ps(ps)

  grid_i = gphys.grid_copy.axis(0)
  grid_j = gphys.grid_copy.axis(2)
  im = gphys.grid_copy.axis(0).pos.val.size
  jm = gphys.grid_copy.axis(2).pos.val.size
  km = gphys.grid_copy.axis(1).pos.val.size
  grid_k =
    Axis.new().
    set_pos(VArray.new(plev).rename("plev").
              put_att("units","hPa"))

  wkgp1 = NArray.sfloat(km).fill(0)
  wkgarr1 = NArray.sfloat(km).fill(0)
  gphys_out = NArray.sfloat(im,plev.size,jm).fill(0)

  jm.times {|j|
    im.times {|i|

      if plev_org[i,0,j] <= plev_org[i,1,j] 
        km.times {|k|
          wkgp1[ k ]   = plev_org[ i, k, j ]
          wkgarr1[ k ] = gphys[ i, k, j ].data.val.to_na
        }
      else
        km.times {|k|
          wkgp1[ k ]   = plev_org[ i, km-(k-1), j ]
          wkgarr1[ k ] = gphys[ i, km-(k-1), j ].data.val.to_na
        }
      end 
        
      plev.size.times { |k|
        
        if ( wkgp1[ 0 ] > plev[ k ] ) then
          gphys_out[ i, k, j ] = -999.0
        else          
          l = km
          1.upto(km-1) {|ll|
#            p "#{wkgp1[ ll ]}, #{plev[ k ]}"
            if ( wkgp1[ ll ] >= plev[ k ] ) then
              l = ll
              break
            end 
          }
          
          if( l == km ) then
            gphys_out[ i, k, j] = -999.0
          else
            gphys_out[ i, k, j ] = 
              ( wkgarr1[ l ] - wkgarr1[ l-1 ] )/( log( wkgp1[ l ]/wkgp1[ l-1 ] ) ) * ( log( plev[ k ] / wkgp1[ l-1 ] ) ) + wkgarr1[ l-1 ]

#             p gphys_out[ i, k, j ]
          end 
       end 
        }

#      p "#{i}, #{j} #{gphys_out[i,j,true].to_a}"
#      p "#{i}, #{j}"
}
}

  grid = Grid.new(grid_i,grid_k,grid_j)
  gphys_out = GPhys.new(grid, VArray.new(gphys_out).rename(gphys.name).
              put_att("units",gphys.data.get_att("units")).
              put_att("ape_name",gphys.data.get_att("ape_name")) )

return gphys_out

end


def sigma_ps(ps)

   a = NArray[ 0.0,
    20.000000,
    38.425343,
    63.647804,
    95.636963,
   134.483307,
   180.584351,
   234.779053,
   298.495789,
   373.971924,
   464.618134,
   575.651001,
   713.218079,
   883.660522,
  1094.834717,
  1356.474609,
  1680.640259,
  2082.273926,
  2579.888672,
  3196.421631,
  3960.291504,
  4906.708496,
  6018.019531,
  7306.631348,
  8765.053711,
 10376.126953,
 12077.446289,
 13775.325195,
 15379.805664,
 16819.474609,
 18045.183594,
 19027.695313,
 19755.109375,
 20222.205078,
 20429.863281,
 20384.480469,
 20097.402344,
 19584.330078,
 18864.750000,
 17961.357422,
 16899.468750,
 15706.447266,
 14411.124023,
 13043.218750,
 11632.758789,
 10209.500977,
  8802.356445,
  7438.803223,
  6144.314941,
  4941.778320,
  3850.913330,
  2887.696533,
  2063.779785,
  1385.912598,
   855.361755,
   467.333588,
   210.393890,
    65.889244,
     7.367743,
     0.000000,
     0.000000] 

  b = NArray[ 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000758235,
 0.0004613950,
 0.0018151561,
 0.0050811190,
 0.0111429105,
 0.0206778757,
 0.0341211632,
 0.0516904071,
 0.0735338330,
 0.0996746942,
 0.1300225109,
 0.1643843204,
 0.2024759352,
 0.2439331412,
 0.2883229554,
 0.3351548910,
 0.3838921487,
 0.4339629412,
 0.4847715795,
 0.5357099175,
 0.5861684084,
 0.6355474591,
 0.6832686067,
 0.7287858129,
 0.7715966105,
 0.8112534285,
 0.8473749161,
 0.8796569109,
 0.9078838825,
 0.9319403172,
    0.9518215060,
    0.9676452279,
    0.9796627164,
    0.9882701039,
    0.9940194488,
    0.9976301193,
    1.0000000000 ]

  data = NArray.sfloat(ps[true,0].data.val.size, 47 ,ps[0,true].data.val.size).fill(0)

  47.times{ |num|
    numnum = num + 14
  data[true,num,true] = a[numnum] + b[numnum]*ps.data.val.to_na
  }

  return data

end