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


# ["T39L48_ias","T39L48_ksc","T39L48_kuo","T39L48_mca","T39L48_non",
#    "T39L24_eml","T39L96_eml","T39L24_non","T39L96_non",
#    "T79L48_eml"].each{ |resol|
#
#    $sstid = "control"
#    $resol = resol
#    mknc_stspctrum
#  }

#  gphys = __gr_tr_lat0("Q")
#  gphys = __gr_tr_lat0("T")
#  gphys = __gr_tr_lat0("U")
#  gphys = __gr_tr_lat0("OMG")

  setopt
  if ARGV.index("-marge")
    mknc_file_marge
  else
    mknc_spctfilter
  end

}

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

  num = 0

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

  $rezol = 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 "$rezol = #{$rezol}, $var = #{$var}, #{Time.now}\n"

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| 

      unless name == "x"  || name == "z"  || name == "time" 
        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
  
  file = NetCDF.create("#{$rezol}_#{$sstid}_#{$filter}filt_#{$var.downcase}.nc")
  puts "#{$rezol}_#{$sstid}_#{$filter}filt_#{$var.downcase}.nc create start"

  $gr2ncfile_path = "/home/yukiko/work/ape/yukiko/data/#{$rezol}_expID01/"
  Dir.foreach($gr2ncfile_path) { |item|
    @data = Ape.new("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/
  }
  
  if $var == "tr_tppn"
    gphys = __gr_tr_tppn
  elsif $var == "tr_tconv"
    gphys = __gr_tr_lat0("pf_t_conv").rename("tr_tconv")
  else
    gphys = __gr_tr_lat0($var)
  end

  if $filter =~ /^k/
    gphys_out = kelvin_filt(gphys)  
  elsif $filter =~ /^a/
    gphys_out = advect_filt(gphys)  
  elsif $filter =~ /^g/
    gphys_out = grav_filt(gphys)  
  end

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

  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 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 = 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 = spct.fft(true).real

#    if $rezol == "T159L48_eml"
#      gphys = gphys[true,-801..-401]
#    else
#      gphys = gphys[true,-401..-1]
#    end
  else

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

      print "lev=#{levnum}, #{Time.now}\n"
      
      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)[levnum,true].data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up)[levnum,true].data.to_a +
          (spct_tmp.cut(num,false,up..999)[levnum,true].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)[levnum,true].data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up)[levnum,true].data.to_a +
          (spct_tmp.cut(num,false,up..999)[levnum,true].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
      
    }

#    if $rezol == "T159L48_eml"
#      gphys = gphys[true,true,-801..-401]
#    else
#      gphys = gphys[true,true,-401..-1]
#    end
  end

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

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

def kelvin_filt(gphys)  

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

  $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)

  if gphys.rank == 2

    spct = spct*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,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

#    if $rezol == "T159L48_eml"
#      gphys = gphys[true,-801..-401]
#    else
#      gphys = gphys[true,-401..-1]
#    end

  else

    levsize = spct.grid_copy.axis(1).pos.val.size
    levsize.times{ |levnum|
      
      print "lev=#{levnum}, #{Time.now}\n"
      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)[levnum,true].data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up)[levnum,true].data.to_a +
          (spct_tmp.cut(num,false,up..999)[levnum,true].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)[levnum,true].data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up)[levnum,true].data.to_a +
          (spct_tmp.cut(num,false,up..999)[levnum,true].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
      
    }

#    if $rezol == "T159L48_eml"
#      gphys = gphys[true,true,-801..-401]
#    else
#      gphys = gphys[true,true,-401..-1]
#    end

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

end


def grav_filt(gphys)  

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

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

  if gphys.rank == 2

    spct = spct*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,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

#    if $rezol == "T159L48_eml"
#      gphys = gphys[true,-801..-401]
#    else
#      gphys = gphys[true,-401..-1]
#    end

  else
    
    levsize = spct.grid_copy.axis(1).pos.val.size
    levsize.times{ |levnum|
      
      print "lev=#{levnum}, #{Time.now}\n"

      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)[levnum,true].data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up)[levnum,true].data.to_a +
          (spct_tmp.cut(num,false,up..999)[levnum,true].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)[levnum,true].data*0.0).to_a +
          spct_tmp.cut(num,false,bo..up)[levnum,true].data.to_a +
          (spct_tmp.cut(num,false,up..999)[levnum,true].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
      
    }

#    if $rezol == "T159L48_eml"
#      gphys = gphys[true,true,-801..-401]
#    else
#      gphys = gphys[true,true,-401..-1]
#    end

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

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 __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 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