#!/usr/bin/env ruby

# -----------------------------------------------
# Ape NetCDF お絵描き応用クラス

class Ape_mkfig

  # class method ------------------------------------------

  def Ape_mkfig.help
    print  <<EOF

  ==========================================================
   class Ape_mkfig のメソッド (ape-view.rb) 
  ---------------------------------------------------------
   * Ape_mkfig.new(num)
   * mkfig_plot(gphys)
   * mkfig_plot_all
   * nc_varlist(file)
   * nc_gt, nc_mf, nc_ml, nc_sh, nc_sh_zonal, nc_tr .....
  ==========================================================

EOF

  end

  # method -----------------------------------------------

  def initialize(num)
    @fig_num = num
  end

  def help
    Ape_mkfig.help
  end

  def nc_varlist(file)
    @data = Ape.new(file)
    @data.netcdf_open.var_names.each { |item|
      unless item =~ /lon|lat|plev|time/ 
	gphys = @data.gphys_open(item)
	var_name = gphys.data.get_att("ape_name").gsub("_", " ")
	print "#{gphys.name}: #{var_name}\n"
      end
    }
  end

  def mkfig_plot(gphys, gphys_u=nil, gphys_v=nil)

    if @fig_num == 3 then
      if gphys_u == nil || gphys_v == nil then
	@data.mkdump(gphys)
      else
        @data.mkdump(gphys, gphys_u, gphys_v)
      end
      
    elsif @fig_num == 2
      
      if gphys_u == nil || gphys_v == nil then
        @data.mkps(gphys)
      else
        @data.mkps(gphys, gphys_u, gphys_v)
      end

    else
      
      if gphys_u == nil || gphys_v == nil then
        @data.plot(gphys)
      else
	@data.plot(gphys, gphys_u, gphys_v)
      end

    end

  end

  def nc_plot_all
    $file_name = nil ;  $pre_file = nil

#    nc_sh
#    nc_tr
#    nc_ml
#    nc_gt
#    nc_mf
#    nc_sh_zonal


=begin
#    nc_sh_anm
    nc_sh_anm
    nc_gt_anm
    nc_ml_anm
    nc_mf_anm
    nc_sh_zonal_anm
    nc_sh
    nc_gt
    nc_ml
    nc_mf
    nc_sh_zonal
    nc_tr
#    nc_tr
=end
  end


  def nc_sst
    $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/sstncfiles/"

    rezol = ["T39", "T79", "T159", "T319"]
    rezol.each{ |item|
      $rezol = item
      $file_label = $rezol
      @data = Ape.new("#{$ncfile_path}#{$rezol}_SST.nc")
      id = ["control","flat","peaked","Qobs","control-5N"]
      gphys_ary = []
      id.each { |expID| 
	$expID = expID
	gphys = @data.gphys_open("sst_#{expID}")
	gphys_ary.push(gphys.
		       set_att("ape_name","sea surcface temperature").
		       set_att("line_name","sst_#{$expID}").mean(0).
		       rename("sst_zonal")
		       )
      }
      mkfig_plot(gphys_ary)

      gphys_con = @data.gphys_open("sst_control") 
      id = ["1keq","3keq","3kw1"]
      id.each { |expID| 
	$expID = expID
	$file_label = "#{$rezol}_#{$expID}"
	gphys = @data.gphys_open("sst_#{expID}") 
	mkfig_plot( (gphys - gphys_con).
		   set_att("ape_name","sea surcface temperature").
		   set_att("line_name","sst_#{$expID}").lon_lotate.
		   set_lost_axes("(diff) from (mean) zonal of control"). 
		   rename("sst_anm")
		   )
      }
      
    }

    id = ["control","flat","peaked","Qobs","control-5N"]
    id.each { |expID| 
      $expID = expID
      $file_label = $expID
      gphys_ary = []
      rezol = ["T39", "T79", "T159", "T319"]
      rezol.each{ |item|
	$rezol = item
	@data = Ape.new("#{$ncfile_path}#{$rezol}_SST.nc")
	gphys = @data.gphys_open("sst_#{expID}")
	gphys_ary.push(gphys.
		       set_att("ape_name","sea surcface temperature").
		       set_att("line_name",$rezol).mean(0).
		       rename("sst_zonal")
		       )
      }
      mkfig_plot(gphys_ary)
    }
  end

  # AGUforAPE-03a_GT_control.nc
  def nc_gt
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item| 
      unless item =~ "time"  || item =~ "gt_slwu" # gt_slwu は一定値
	gphys = @data.gphys_open(item)
	mkfig_plot(gphys)
      end
    }

    # total_precipitation_flux
    g1 = @data.gphys_open("gt_cppn")
    g2 = @data.gphys_open("gt_dppn")
    gphys = (g1 + g2).rename("gt_tppn").
      set_att("ape_name", "total_precipitation_flux")
    mkfig_plot(gphys)

  end

  # #{$groupid}_GT_control.nc
  def nc_ga_gt
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
    print "#{$rezol}_#{$expID}\n"
    @data.netcdf_open.var_names.each { |item| 
      unless item =~ "time"    # gt_slwu は一定値
	gphys = @data.gphys_open(item)
	print "#{gphys.name}: #{gphys.mean(0)}\n"
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("gt_cppn")
    g2 = @data.gphys_open("gt_dppn")
    gphys = (g1 + g2).rename("gt_tppn").
      set_att("ape_name", "total_precipitation_flux")
    print "#{gphys.name}: #{gphys.mean(0)}\n"
  end

  # #{$groupid}_ML_control.nc
  def nc_ml
    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "lon" || item =~ "plev"  
	gphys = @data.gphys_open(item)
	mkfig_plot(gphys.mean(0))
      end
    }
    # mass_stream_function
    gphys = @data.gphys_open("ml_v")
    gphys = gphys.mean(0).strm_function("ml_psi").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1") 
    mkfig_plot(gphys)
  end

  # #{$groupid}_TR_control.nc
  def nc_tr
    @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |name|
      if name =~ /lw_toa|mslp|500|250/
	if name =~ /lw_toa|mslp/
	  lost_axis = []
	elsif name =~ /250/
	  lost_axis = ["plev=25000 Pa"]
	elsif name =~ /500/
	  lost_axis = ["plev=50000 Pa"]
	end

	# x-t ダイヤグラム
	gphys = @data.gphys_open(name)
	gphys =  gphys.cut(true,0,true)[true,-400..-1]
#	mkfig_plot((gphys - gphys.mean(0)).
#		   add_lost_axes(lost_axis).
#		   add_lost_axes("(diff) from (mean) zonal").lon_lotate)

	# 時空間スペクトル
	gphys = gphys.stspct_fft("#{name}_spct").
	  set_att("ape_name",
		  "#{gphys.data.get_att("ape_name")}_S-T_power_spectrum").
	  set_att("units","(#{gphys.data.get_att("units")} day-1)2").
	  add_lost_axes(gphys.lost_axes).
	  add_lost_axes(lost_axis)
	dim1 = gphys.coord(0).shape.to_s.to_i
	dim2 = gphys.coord(1).shape.to_s.to_i
#	mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..50])

      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
#	gphys = gphys.cut(true,0,true)[true,-400..-1]
	gphys = gphys.cut(true,0,true)[true,-401..-1]

#=begin
	grid_0 = gphys.grid_copy.axis(0)
	grid_1 =
	  Axis.new().
	  set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time").
 put_att("units","days"))
        grid = Grid.new(grid_0,grid_1)
        gphys = GPhys.new(grid, gphys.data )
#=end

#	mkfig_plot(gphys.lon_lotate)
	mkfig_plot(gphys.lon_lotate.rename("tr_tppn_mono").set_lost_axes(""))

	# 時空間スペクトル
	gphys = gphys.stspct_fft("tr_tppn_spct").
	  set_att("ape_name",
		  "#{gphys.data.get_att("ape_name")}_S-T_power_spectrum").
	  set_att("units","(#{gphys.data.get_att("units")} day-1)2").
	  add_lost_axes(gphys.lost_axes)
	dim1 = gphys.coord(0).shape.to_s.to_i
	dim2 = gphys.coord(1).shape.to_s.to_i
#	mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..50])
      end
    }

  end

  #  #{$groupid}_SH_control.nc
  def nc_sh
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys = @data.gphys_open(item)
	gphys = gphys.lon_lotate
	mkfig_plot(gphys)
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2).rename("sh_tppn").
      set_att("ape_name", "total_precipitation_flux")
    gphys = gphys.lon_lotate
    mkfig_plot(gphys)
  end

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys = @data.gphys_open(item)
	gphys = gphys.rename("#{gphys.name}_zonal")
	mkfig_plot(gphys.mean(0))
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2).rename("sh_tppn_zonal").
      set_att("ape_name", "total_precipitation_flux")
    mkfig_plot(gphys.mean(0))
  end
  
  #  #{$groupid}_MF_control.nc
  def nc_mf
    @data = Ape.new("#{$ncfile_path}#{$groupid}_MF_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "plev" 
	gphys = @data.gphys_open(item)
	mkfig_plot(gphys)
      end
    }
  end

  # anomary -----------------------------------------------------

  def nc_gt_anm(zonalcont="control")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{zonalcont}.nc")
    @data.netcdf_open.var_names.each { |item| 
      unless item =~ "time"  || item =~ "gt_slwu" # gt_slwu は一定値
	gphys_con = @data_con.gphys_open(item).mean(0)
	lost_axis = ["(diff) from (mean) time of #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys - gphys_con)
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("gt_cppn")
    g2 = @data.gphys_open("gt_dppn")
    gphys = (g1 + g2)
    g1 = @data_con.gphys_open("gt_cppn")
    g2 = @data_con.gphys_open("gt_dppn")
    gphys_con = (g1 + g2).mean(0)
    lost_axis = ["(diff) from (mean) time of #{zonalcont}"]
    gphys = (gphys - gphys_con).rename("gt_tppn_anm").
      set_att("ape_name", "total_precipitation_flux").
      add_lost_axes(lost_axis)
    mkfig_plot(gphys)
  end

  def nc_ml_anm(zonalcont="control")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "lon" || item =~ "plev"  
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys.mean(0) - gphys_con.mean(0))
      end
    }
    # mass_stream_function
    lost_axis = ["(diff) from (mean) #{zonalcont}"]
    gphys = @data.gphys_open("ml_v")
    gphys = gphys.add_lost_axes(lost_axis).mean(0).strm_function("ml_psi_anm").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1")
    gphys_con = @data_con.gphys_open("ml_v")
    gphys_con = gphys_con.mean(0).strm_function("ml_psi_anm").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1") 
    mkfig_plot(gphys-gphys_con)
  end

  def nc_sh_anm(zonalcont="control")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{zonalcont}.nc")

    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) zonal of #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = (gphys - gphys_con.mean(0))
	gphys = gphys.lon_lotate.rename("#{gphys.name}_anm")
	mkfig_plot(gphys)
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2)
    g1 = @data_con.gphys_open("sh_cppn")
    g2 = @data_con.gphys_open("sh_dppn")
    gphys_con = (g1 + g2)
    gphys = (gphys - gphys_con.mean(0))
    lost_axis = ["(diff) from (mean) zonal of #{zonalcont}"]
    gphys = gphys.lon_lotate.rename("sh_tppn_anm").
      set_att("ape_name", "total_precipitation_flux").
      add_lost_axes(lost_axis)
#    mkfig_plot(gphys)

    # 赤道上 降水量 (tppn) 
    gphys = gphys.set_lost_axes("")
    gphys =  gphys.cut(true,0).rename("sh_tppn_lat0_anm"). 
      add_lost_axes(lost_axis)
    mkfig_plot(gphys)
    # 赤道上 蒸発量 (slh) 
    gphys = @data.gphys_open("sh_slh")
    gphys_con = @data_con.gphys_open("sh_slh")
    gphys = (gphys - gphys_con.mean(0)).lon_lotate
    gphys = gphys.cut(true,0).rename("sh_slh_lat0_anm").
      add_lost_axes(lost_axis)
    mkfig_plot(gphys)
    
  end

  def nc_sh_zonal_anm(zonalcont="control")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{zonalcont}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_zonal_anm")
	mkfig_plot(gphys.mean(0)- gphys_con.mean(0))
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2)
    g1 = @data_con.gphys_open("sh_cppn")
    g2 = @data_con.gphys_open("sh_dppn")
    gphys_con = (g1 + g2)
    gphys = (gphys.mean(0) - gphys_con.mean(0))
    lost_axis = ["(diff) from (mean) #{zonalcont}"]
    gphys = gphys.rename("sh_tppn_zonal_anm").
      set_att("ape_name", "total_precipitation_flux").
      add_lost_axes(lost_axis)
    mkfig_plot(gphys)
  end
  
  def nc_mf_anm(zonalcont="control")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_MF_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_MF_#{zonalcont}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "plev" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys - gphys_con)
      end
    }
  end

  # Hosaka plot --------------------------------------

  def nc_sh_psuv_anm(zonalcont="control")

    lost_axis = ["plev=100000 Pa",
      "(diff) from (mean) #{zonalcont}"]
   
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{zonalcont}.nc")
    gphys_ps = @data.gphys_open("sh_ps")
    gphys_ps_con = @data_con.gphys_open("sh_ps")

    gphys_ps = ( gphys_ps - gphys_ps_con.mean(0) ).lon_lotate.
      add_lost_axes(lost_axis).
      rename("sh_psuv_anm").
      set_att("ape_name","ps, (u,v)").
      set_att("units","Pa, m s-1")

    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")
    gphys_u = @data.gphys_open("ml_u").cut(true,true,100000)
    gphys_v = @data.gphys_open("ml_v").cut(true,true,100000)
    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,100000)
    gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,100000)

    gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
    gphys_v = ( gphys_v - gphys_v_con.mean(0) ).lon_lotate


    u_data  = NArray.sfloat( gphys_u.axis(0).pos.val.size ,
			    gphys_u.axis(1).pos.val.size ).fill!(0.0)
    v_data  = NArray.sfloat( gphys_v.axis(0).pos.val.size ,
			    gphys_v.axis(1).pos.val.size ).fill!(0.0)
    
    num_x = gphys_u.axis(0).pos.val.size/40
    num_y = gphys_u.axis(1).pos.val.size/20

    gphys_u.axis(0).pos.val.size.times{ |num|
      if (num % num_x) == 0
	u_data[num,true]  =  gphys_u.data.val[num,true]
	v_data[num,true]  =  gphys_v.data.val[num,true]
      elsif   
	u_data[num,true]  = 0.0
	v_data[num,true]  = 0.0
      end
    }  
    
    gphys_u.axis(1).pos.val.size.times{ |num|
      if ((num+2) % num_y) == 0
	u_data[true,num]  =  u_data[true,num]
	v_data[true,num]  =  v_data[true,num]
      elsif
	u_data[true,num]  = 0.0
	v_data[true,num]  = 0.0
      end
    }

    mkfig_plot(gphys_ps,u_data,v_data)
    
  end
  
  def nc_sh_phiuv250_anm(zonalcont="control")

    lost_axis = ["(diff) from (mean) #{zonalcont}"]

    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")

    gphys_phi = @data.gphys_open("ml_phi").cut(true,true,25000)
    gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,25000)

    gphys_phi = ( gphys_phi - gphys_phi_con.mean(0) ).lon_lotate.
      add_lost_axes(lost_axis).
      rename("sh_phiuv250_anm").
      set_att("ape_name","phi, (u,v)").
      set_att("units","m, m s-1")

    gphys_u = @data.gphys_open("ml_u").cut(true,true,25000)
    gphys_v = @data.gphys_open("ml_v").cut(true,true,25000)
    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,25000)
    gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,25000)

    gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
    gphys_v = ( gphys_v - gphys_v_con.mean(0) ).lon_lotate


    u_data  = NArray.sfloat( gphys_u.axis(0).pos.val.size ,
			    gphys_u.axis(1).pos.val.size ).fill!(0.0)
    v_data  = NArray.sfloat( gphys_v.axis(0).pos.val.size ,
			    gphys_v.axis(1).pos.val.size ).fill!(0.0)
    
    num_x = gphys_u.axis(0).pos.val.size/40
    num_y = gphys_u.axis(1).pos.val.size/20

    gphys_u.axis(0).pos.val.size.times{ |num|
      if (num % num_x) == 0
	u_data[num,true]  =  gphys_u.data.val[num,true]
	v_data[num,true]  =  gphys_v.data.val[num,true]
      elsif   
	u_data[num,true]  = 0.0
	v_data[num,true]  = 0.0
      end
    }  
    
    gphys_u.axis(1).pos.val.size.times{ |num|
      if ((num+2) % num_y) == 0
	u_data[true,num]  =  u_data[true,num]
	v_data[true,num]  =  v_data[true,num]
      elsif
	u_data[true,num]  = 0.0
	v_data[true,num]  = 0.0
      end
    }

    mkfig_plot(gphys_phi,u_data,v_data)
    
  end
  
  def nc_sh_phiuv850_anm(zonalcont="control")

    lost_axis = ["(diff) from (mean) #{zonalcont}"]

    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")

    gphys_phi = @data.gphys_open("ml_phi").cut(true,true,85000)
    gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,85000)

    gphys_phi = ( gphys_phi - gphys_phi_con.mean(0) ).lon_lotate.
      add_lost_axes(lost_axis).
      rename("sh_phiuv850_anm").
      set_att("ape_name","phi, (u,v)").
      set_att("units","m, m s-1")

    gphys_u = @data.gphys_open("ml_u").cut(true,true,85000)
    gphys_v = @data.gphys_open("ml_v").cut(true,true,85000)
    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,85000)
    gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,85000)

    gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
    gphys_v = ( gphys_v - gphys_v_con.mean(0) ).lon_lotate


    u_data  = NArray.sfloat( gphys_u.axis(0).pos.val.size ,
			    gphys_u.axis(1).pos.val.size ).fill!(0.0)
    v_data  = NArray.sfloat( gphys_v.axis(0).pos.val.size ,
			    gphys_v.axis(1).pos.val.size ).fill!(0.0)
    
    num_x = gphys_u.axis(0).pos.val.size/40
    num_y = gphys_u.axis(1).pos.val.size/20

    gphys_u.axis(0).pos.val.size.times{ |num|
      if (num % num_x) == 0
	u_data[num,true]  =  gphys_u.data.val[num,true]
	v_data[num,true]  =  gphys_v.data.val[num,true]
      elsif   
	u_data[num,true]  = 0.0
	v_data[num,true]  = 0.0
      end
    }  
    
    gphys_u.axis(1).pos.val.size.times{ |num|
      if ((num+2) % num_y) == 0
	u_data[true,num]  =  u_data[true,num]
	v_data[true,num]  =  v_data[true,num]
      elsif
	u_data[true,num]  = 0.0
	v_data[true,num]  = 0.0
      end
    }

    mkfig_plot(gphys_phi,u_data,v_data)
    
  end
  
  def nc_sh_tuv500_anm(zonalcont="control")

    lost_axis = ["(diff) from (mean) #{zonalcont}"]

    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")

    gphys_t = @data.gphys_open("ml_t").cut(true,true,50000)
    gphys_t_con = @data_con.gphys_open("ml_t").cut(true,true,50000)

    gphys_t = ( gphys_t - gphys_t_con.mean(0) ).lon_lotate.
      add_lost_axes(lost_axis).
      rename("sh_tuv500_anm").
      set_att("ape_name","T, (u,v)").
      set_att("units","K, m s-1")

    gphys_u = @data.gphys_open("ml_u").cut(true,true,50000)
    gphys_v = @data.gphys_open("ml_v").cut(true,true,50000)
    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,50000)
    gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,50000)

    gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
    gphys_v = ( gphys_v - gphys_v_con.mean(0) ).lon_lotate


    u_data  = NArray.sfloat( gphys_u.axis(0).pos.val.size ,
			    gphys_u.axis(1).pos.val.size ).fill!(0.0)
    v_data  = NArray.sfloat( gphys_v.axis(0).pos.val.size ,
			    gphys_v.axis(1).pos.val.size ).fill!(0.0)
    
    num_x = gphys_u.axis(0).pos.val.size/40
    num_y = gphys_u.axis(1).pos.val.size/20

    gphys_u.axis(0).pos.val.size.times{ |num|
      if (num % num_x) == 0
	u_data[num,true]  =  gphys_u.data.val[num,true]
	v_data[num,true]  =  gphys_v.data.val[num,true]
      elsif   
	u_data[num,true]  = 0.0
	v_data[num,true]  = 0.0
      end
    }  
    
    gphys_u.axis(1).pos.val.size.times{ |num|
      if ((num+2) % num_y) == 0
	u_data[true,num]  =  u_data[true,num]
	v_data[true,num]  =  v_data[true,num]
      elsif
	u_data[true,num]  = 0.0
	v_data[true,num]  = 0.0
      end
    }

    mkfig_plot(gphys_t,u_data,v_data)
    
  end

  
  def nc_ml_tuom_anm(zonalcont="control")

    lost_axis = ["(diff) from (mean) #{zonalcont}"]

    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")

    gphys_t = @data.gphys_open("ml_t").cut(true,0,true)
    gphys_t_con = @data_con.gphys_open("ml_t").cut(true,0,true)

    gphys_t = ( gphys_t - gphys_t_con.mean(0) ).lon_lotate.
      add_lost_axes(lost_axis).
      rename("ml_tuom_anm").
      set_att("ape_name","T, (u,-omg)").
      set_att("units","K, m s-1")

    gphys_u = @data.gphys_open("ml_u").cut(true,0,true)
    gphys_om = @data.gphys_open("ml_v").cut(true,0,true)
    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,0,true)
    gphys_om_con = @data_con.gphys_open("ml_v").cut(true,0,true)

    gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
    gphys_om = ( gphys_om - gphys_om_con.mean(0) ).lon_lotate


    u_data  = NArray.sfloat( gphys_u.axis(0).pos.val.size ,
			    gphys_u.axis(1).pos.val.size ).fill!(0.0)
    om_data  = NArray.sfloat( gphys_om.axis(0).pos.val.size ,
			    gphys_om.axis(1).pos.val.size ).fill!(0.0)
    
    num_x = gphys_u.axis(0).pos.val.size/40

    gphys_u.axis(0).pos.val.size.times{ |num|
      if (num % num_x) == 0
	u_data[num,true]  =  gphys_u.data.val[num,true]
	om_data[num,true]  =  gphys_om.data.val[num,true]
      elsif   
	u_data[num,true]  = 0.0
	om_data[num,true]  = 0.0
      end
    }  
    
    mkfig_plot(gphys_t,u_data,-om_data)
    
  end

  
  # line 重ね描き -----------------------------------------------------

  # #{$groupid}_GT_control.nc
  def nc_gt_comparing_expid
    id = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
#    id = ["control","flat","peaked","Qobs","control-5N"]
    $file_label = "#{$rezol}"
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item =~ "time"  || item =~ "gt_slwu" # gt_slwu は一定値
	id.each{ |expID|
	  $expID = expID
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
	  gphys_ary.push(@data.gphys_open(item).
			 set_att("line_name",$expID) )
	}
	mkfig_plot(gphys_ary)
      end
    }

    # total_precipitation_flux
    gphys_ary = []
    id.each{ |expID|
      $expID = expID
      @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
      g1 = @data.gphys_open("gt_cppn")
      g2 = @data.gphys_open("gt_dppn")
      gphys = (g1 + g2).rename("gt_tppn").
	set_att("ape_name", "total_precipitation_flux")
      gphys_ary.push(gphys.
		     set_att("line_name",$expID) )
    }
    mkfig_plot(gphys_ary)
  end
 
  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal_comparing_expid
    id = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
#    id = ["control","flat","peaked","Qobs","control-5N"]
    $file_label = "#{$rezol}"
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item =~ "lon" || item =~ "lat" 
	id.each{ |expID|
	  $expID = expID
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
	  gphys = @data.gphys_open(item).mean(0)
  	  gphys_ary.push(gphys.rename("#{gphys.name}_zonal").
			 set_att("line_name",$expID) )
	}
	mkfig_plot(gphys_ary)
      end
    }
    # total_precipitation_flux
    gphys_ary = []
    id.each{ |expID|
      $expID = expID
      @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
      g1 = @data.gphys_open("sh_cppn")
      g2 = @data.gphys_open("sh_dppn")
      gphys = (g1 + g2).rename("sh_tppn_zonal").
	set_att("ape_name", "total_precipitation_flux")
      gphys_ary.push(gphys.mean(0).
		     set_att("line_name",$expID) )
    }
    mkfig_plot(gphys_ary)
  end

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_anm_comparing_expid
    id_con = ["control","flat","Qobs","H1998con"]
    id     = ["3keq","flat3keq","Qobs3keq","H1998pa"]
    $file_label = "#{$rezol}"
    
    gphys_ary_tppn = []
    gphys_ary_slh = []
    id.size.times{ |num|
      @data_con = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{id_con[num]}.nc")
      @data     = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{id[num]}.nc")

      # 赤道上 降水量 (tppn) 
      g1 = @data_con.gphys_open("sh_cppn")
      g2 = @data_con.gphys_open("sh_dppn")
      gphys_con = (g1 + g2).rename("sh_tppn_zonal")
      g1 = @data.gphys_open("sh_cppn")
      g2 = @data.gphys_open("sh_dppn")
      gphys= (g1 + g2)
      gphys = (gphys - gphys_con.mean(0)).lon_lotate
      gphys =  gphys.cut(true,0).rename("sh_tppn_lat0_anm"). 
	set_att("ape_name", "total_precipitation_flux"). 
	set_att("line_name","#{id_con[num]}-#{id[num]}")
      gphys_ary_tppn.push(gphys)
      
      # 赤道上 蒸発量 (slh) 
      gphys = @data.gphys_open("sh_slh")
      gphys_con = @data_con.gphys_open("sh_slh")
      gphys = (gphys - gphys_con.mean(0)).lon_lotate
      gphys = gphys.cut(true,0).rename("sh_slh_lat0_anm").
	set_att("line_name","#{id_con[num]}-#{id[num]}")
      gphys_ary_slh.push(gphys)

    }
    mkfig_plot(gphys_ary_tppn)
    mkfig_plot(gphys_ary_slh)

  end


  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal_comparing_rezol
    $file_label = "#{$cumulus}_#{$expID}"
#    rezol = ["T39L48_eml", "T39L24_eml", "T79L48_eml", "T159L48_eml"]
    rezol = ["T39L48_eml", "T39L24_eml", "T79L48_eml"]
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item =~ "lon" || item =~ "lat" 
	rezol.each{ |rezol_item|
	  $rezol = rezol_item
	  $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
	  gphys = @data.gphys_open(item).mean(0)
  	  gphys_ary.push(gphys.rename("#{gphys.name}_zonal").
			 set_att("line_name",$rezol) )
	}
	mkfig_plot(gphys_ary)
      end
    }
    # total_precipitation_flux
    gphys_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
      g1 = @data.gphys_open("sh_cppn")
      g2 = @data.gphys_open("sh_dppn")
      gphys = (g1 + g2).rename("sh_tppn_zonal").
	set_att("ape_name", "total_precipitation_flux")
      gphys_ary.push(gphys.mean(0).
		     set_att("line_name",$rezol) )
    }
    mkfig_plot(gphys_ary)
  end

  # #{$groupid}_GT_control.nc
  def nc_gt_comparing_rezol
    $file_label = "#{$cumulus}_#{$expID}"
#    rezol = ["T39L48_eml", "T39L24_eml", "T79L48_eml", "T159L48_eml"]
    rezol = ["T39L48_eml", "T39L24_eml", "T79L48_eml"]
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item =~ "time"  || item =~ "gt_slwu" # gt_slwu は一定値
	rezol.each{ |rezol_item|
	  $rezol = rezol_item
	  $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
	  gphys_ary.push(@data.gphys_open(item).
			 set_att("line_name",$rezol) )
	  puts $rezol
	}
	mkfig_plot(gphys_ary)
      end
    }

    # total_precipitation_flux
    gphys_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
      g1 = @data.gphys_open("gt_cppn")
      g2 = @data.gphys_open("gt_dppn")
      gphys = (g1 + g2).rename("gt_tppn").
	set_att("ape_name", "total_precipitation_flux")
      gphys_ary.push(gphys.
		     set_att("line_name",$rezol) )
    }
    mkfig_plot(gphys_ary)
  end



  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_ml_zonal_comparing_rezol
    $file_label = "#{$cumulus}_#{$expID}"
    rezol = ["T39L48_eml", "T39L24_eml","T39L96_eml", "T79L48_eml", "T159L48_eml"]

    # U の赤道鉛直断面
    gphys_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      gphys = @data.gphys_open("ml_u").cut(true,0,true).mean(0)
      gphys_ary.push(gphys.rename("#{gphys.name}_lat0_zonal").
		     set_att("line_name",$rezol) )
    }
    mkfig_plot(gphys_ary)
  end


  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_anm_comparing_rezol

    rezol = ["T39L48_eml", "T39L24_eml","T39L96_eml", "T79L48_eml", "T159L48_eml"]

    expID_con_hash =  {
      "3keq" => "control", 
      "flat3keq" => "flat",
      "Qobs3keq" => "Qobs",
      "H1998pa" => "H1998con"
    }
    $file_label = "#{$cumulus}_#{$expID}"
    zonalcont = expID_con_hash[$expID]
    lost_axis = ["(diff) from (mean) zonal of #{zonalcont}"]

    gphys_ary_tppn = []
    gphys_ary_slh = []

    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      @data_con = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{zonalcont}.nc")
      @data     = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")

      # 赤道上 降水量 (tppn) 
      g1 = @data_con.gphys_open("sh_cppn")
      g2 = @data_con.gphys_open("sh_dppn")
      gphys_con = (g1 + g2).rename("sh_tppn_zonal")
      g1 = @data.gphys_open("sh_cppn")
      g2 = @data.gphys_open("sh_dppn")
      gphys= (g1 + g2)
      gphys = (gphys - gphys_con.mean(0)).lon_lotate
      gphys =  gphys.cut(true,0).rename("sh_tppn_lat0_anm"). 
	set_att("ape_name", "total_precipitation_flux"). 
	set_att("line_name","#{$rezol}").
	add_lost_axes(lost_axis) 

      gphys_ary_tppn.push(gphys)
      
      # 赤道上 蒸発量 (slh) 
      gphys = @data.gphys_open("sh_slh")
      gphys_con = @data_con.gphys_open("sh_slh")
      gphys = (gphys - gphys_con.mean(0)).lon_lotate
      gphys = gphys.cut(true,0).rename("sh_slh_lat0_anm").
	set_att("line_name","#{$rezol}").
	add_lost_axes(lost_axis) 
      gphys_ary_slh.push(gphys)

    }
    mkfig_plot(gphys_ary_tppn)
    mkfig_plot(gphys_ary_slh)

  end


  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal_comparing_cumulus
    $file_label = "T39L48_#{$expID}"
#    rezol = ["T39L48_eml", "T39L48_ias", "T39L48_kuo","T39L48_mca","T39L48_non"]
    rezol = ["T39L48_eml", "T39L48_ias", "T39L48_ksc", "T39L48_kuo","T39L48_mca","T39L48_non"]
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item =~ "lon" || item =~ "lat" 
	rezol.each{ |rezol_item|
	  $rezol = rezol_item
	  $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
	  $ncfile_path = "/mnt/yukiko/netCDF/#{$rezol}/"
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
	  gphys = @data.gphys_open(item).mean(0)
  	  gphys_ary.push(gphys.rename("#{gphys.name}_zonal").
			 set_att("line_name",$rezol) )
	}
	mkfig_plot(gphys_ary)
      end
    }
    # total_precipitation_flux
    gphys_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      $ncfile_path = "/mnt/yukiko/netCDF/#{$rezol}/"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
      g1 = @data.gphys_open("sh_cppn")
      g2 = @data.gphys_open("sh_dppn")
      gphys = (g1 + g2).rename("sh_tppn_zonal").
	set_att("ape_name", "total_precipitation_flux")
      gphys_ary.push(gphys.mean(0).
		     set_att("line_name",$rezol) )
    }
    mkfig_plot(gphys_ary)
  end

  # #{$groupid}_GT_control.nc
  def nc_gt_comparing_cumulus
    $file_label = "T39L48_#{$expID}"
    rezol = ["T39L48_eml", "T39L48_ias", "T39L48_ksc", "T39L48_kuo","T39L48_mca","T39L48_non"]
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item =~ "time"  || item =~ "gt_slwu" # gt_slwu は一定値
	rezol.each{ |rezol_item|
	  $rezol = rezol_item
	  $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
	  $ncfile_path = "/mnt/yukiko/netCDF/#{$rezol}/"
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
	  gphys_ary.push(@data.gphys_open(item).
			 set_att("line_name",$rezol) )
	  puts $rezol
	}
	mkfig_plot(gphys_ary)
      end
    }

    # total_precipitation_flux
    gphys_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      $ncfile_path = "/mnt/yukiko/netCDF/#{$rezol}/"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
      g1 = @data.gphys_open("gt_cppn")
      g2 = @data.gphys_open("gt_dppn")
      gphys = (g1 + g2).rename("gt_tppn").
	set_att("ape_name", "total_precipitation_flux")
      gphys_ary.push(gphys.
		     set_att("line_name",$rezol) )
    }
    mkfig_plot(gphys_ary)
  end


  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_anm_comparing_cumulus

    rezol = ["T39L48_eml", "T39L48_ias", "T39L48_ksc", "T39L48_kuo","T39L48_mca","T39L48_non"]

    expID_con_hash =  {
      "3keq" => "control", 
      "flat3keq" => "flat",
      "Qobs3keq" => "Qobs",
      "H1998pa" => "H1998con"
    }

    $file_label = "T39L48_#{$expID}"
    zonalcont = expID_con_hash[$expID]
    lost_axis = ["(diff) from (mean) zonal of #{zonalcont}"]

    gphys_ary_tppn = []
    gphys_ary_slh = []

    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      @data_con = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{zonalcont}.nc")
      @data     = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")

      # 赤道上 降水量 (tppn) 
      g1 = @data_con.gphys_open("sh_cppn")
      g2 = @data_con.gphys_open("sh_dppn")
      gphys_con = (g1 + g2).rename("sh_tppn_zonal")
      g1 = @data.gphys_open("sh_cppn")
      g2 = @data.gphys_open("sh_dppn")
      gphys= (g1 + g2)
      gphys = (gphys - gphys_con.mean(0)).lon_lotate
      gphys =  gphys.cut(true,0).rename("sh_tppn_lat0_anm"). 
	set_att("ape_name", "total_precipitation_flux"). 
	set_att("line_name","#{$rezol}").
	add_lost_axes(lost_axis) 

      gphys_ary_tppn.push(gphys)
      
      # 赤道上 蒸発量 (slh) 
      gphys = @data.gphys_open("sh_slh")
      gphys_con = @data_con.gphys_open("sh_slh")
      gphys = (gphys - gphys_con.mean(0)).lon_lotate
      gphys = gphys.cut(true,0).rename("sh_slh_lat0_anm").
	set_att("line_name","#{$rezol}").
	add_lost_axes(lost_axis) 
      gphys_ary_slh.push(gphys)

    }
    mkfig_plot(gphys_ary_tppn)
    mkfig_plot(gphys_ary_slh)

  end

  # GrADS -----------------------------------------------------

  def gr_tr_all

    file_name = Array.new
    Dir.foreach($gr2ncfile_path) { |item|
      file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/
    }
    @data = Ape.new(file_name)

#    __gr_tr_mslp

    Ape.new(file_name[0]).netcdf_open.var_names.each { |item|

      if item == "tr_tppn" 
	__gr_tr_tppn
      elsif item == "tr_mslp"
#	__gr_tr_mslp
      elsif item == "U"
#	__gr_tr_u250
      elsif item == "V"
#	__gr_tr_v250
      elsif item == "OMG" 
#	__gr_tr_om500
      end
    }

  end

  def gr_comp_tend

    file_name = Array.new
    Dir.foreach($gr2ncfile_path) { |item|
      file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/
    }
    @data = Ape.new(file_name)

    #  tr_tppn
    rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate
    lost_axes = rain.lost_axes.to_s.sub("y=","lat=")
    rain = rain[true,-400..-1].set_lost_axes(lost_axes)

    # t_conv composite
    gphys = @data.gphys_open("pf_t_conv")[true,true,-400..-1].lon_lotate
    comp = gphys.composite(rain, 0.0007)
    comp = comp.      
      set_att("ape_name",
	      "cumulus_heating_composite" ).
      set_att("units", gphys.data.get_att("units") ).
      rename("comp_t_conv").
      set_lost_axes(lost_axes)

    # 軸の long_name 直し (pressure level => sigma)
    z = comp.grid_copy.axis(1).pos.
      set_att("long_name","sigma").
      set_att("units","1")
    z = Axis.new().set_pos(z)
    grid = comp.grid_copy.change_axis(1, z)
    comp = GPhys.new( grid, comp.data ) 

    # 描画
    dim = comp.axis(0).pos.val.size
    mkfig_plot(comp[(dim/3)..(dim*2/3),true])
    comp = comp.cut(0,true).rename("comp_t_conv_zonal")
    mkfig_plot(comp)

  end

  def gr_comp_point

    file_name = Array.new
    Dir.foreach($gr2ncfile_path) { |item|
      file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/
    }
    @data = Ape.new(file_name)

    #  tr_tppn
    rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate
    lost_axes = rain.lost_axes.to_s.sub("y=","lat=")
    rain = rain[true,-400..-1].set_lost_axes(lost_axes)

    # t_conv composite
    gphys = @data.gphys_open("T")[true,true,-400..-1].lon_lotate
    comp = gphys.composite(rain, 0.0007, true)
    comp = comp.      
      set_att("ape_name",
	      "composite_point" ).
      set_att("units", "1" ).
      rename("comp_point").
      set_lost_axes(lost_axes)

    # 描画
    mkfig_plot(comp)

  end

  def gr_comp_tuv

    file_name = Array.new
    Dir.foreach($gr2ncfile_path) { |item|
      file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/
    }
    @data = Ape.new(file_name)

    #  tr_tppn
    rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate
    lost_axes = rain.lost_axes.to_s.sub("y=","lat=")
    rain = rain[true,-400..-1].set_lost_axes(lost_axes)

    # T composite
    gphys = @data.gphys_open("T")[true,true,-400..-1].lon_lotate
    t_comp = gphys.composite(rain, 0.0007)
    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_lost_axes(lost_axes).
      add_lost_axes("(diff) from (mean) zonal").
      rename("comp_tuv")

    # U composite
    gphys = @data.gphys_open("U")[true,true,-400..-1].lon_lotate
    u_comp = gphys.composite(rain, 0.0007)
    u_comp = u_comp - u_comp.mean(0)

    # OMG composite
    gphys = @data.gphys_open("OMG")[true,true,-400..-1].lon_lotate
    om_comp = gphys.composite(rain, 0.0007)    
    om_comp = om_comp - om_comp.mean(0)

    # ベクトルの間引
    u_data  = NArray.sfloat( u_comp.axis(0).pos.val.size , 
			   u_comp.axis(1).pos.val.size ).fill!(0.0)
 
    om_data = NArray.sfloat( om_comp.axis(0).pos.val.size , 
			    om_comp.axis(1).pos.val.size ).fill!(0.0)

    num_l = u_comp.axis(0).pos.val.size/120
    num_z = u_comp.axis(1).pos.val.size/48

    u_comp.axis(0).pos.val.size.times{ |num|
      if (num % num_l) == 0
	u_data[num,true]  =  u_comp.data.val[num,true]
	om_data[num,true] =  om_comp.data.val[num,true]
      elsif
	u_data[num,true]  = 0.0
	om_data[num,true] = 0.0
      end
    }
    
    num_z =1 if num_z < 1

    u_comp.axis(1).pos.val.size.times{ |num|
      if (num % num_z) == 0
	u_data[true,num]  =  u_data[true,num]
	om_data[true,num] =  om_data[true,num]
      elsif
	u_data[true,num]  = 0.0
	om_data[true,num] = 0.0
      end
    }


    # 軸の 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 )

    # 描画
    dim = t_comp.axis(0).pos.val.size
    mkfig_plot(t_comp[(dim/3)..(dim*2/3),true],
	       u_data[(dim/3)..(dim*2/3),true], 
	       -om_data[(dim/3)..(dim*2/3),true])

  end



  # 直接呼び出さない -------------------------------------------
  # ローカルメソッドって __ で書くんだっけ?

  private
  def __gr_tr_tppn
    # tr_tppn
    gphys = @data.gphys_open("tr_tppn").cut(true,0,true)
    lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
#    gphys = gphys[true,-400..-1].set_lost_axes(lost_axes)
    gphys = gphys[true,-401..-1].set_lost_axes(lost_axes)

    grid_0 = gphys.grid_copy.axis(0)
    grid_1 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time").
 put_att("units","days"))
    grid = Grid.new(grid_0,grid_1)
    gphys = GPhys.new(grid, gphys.data )

    mkfig_plot(gphys.lon_lotate.rename("tr_tppn_mono").set_lost_axes(""))

#    mkfig_plot(gphys.lon_lotate)
    # 時空間スペクトル
#    __gr_spct(gphys)
  end

  def __gr_tr_mslp
    # tr_mslp
    gphys = @data.gphys_open("tr_mslp").cut(true,0,true)
    lost_axis = gphys.lost_axes.to_s.sub("y=","lat=")
    gphys = gphys[true,-400..-1].
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate.
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)
  end

  def __gr_tr_u250
    # tr_u250
    gphys = @data.gphys_open("U").cut(true,0.25,true)
    lost_axis = [gphys.data.get_att("lost_axis"),
      gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    gphys = gphys[true,-400..-1].rename("tr_u250").
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate.
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)
  end

  def __gr_tr_v250
    # V
    gphys = @data.gphys_open("V").cut(true,0.25,true)
    lost_axis = [gphys.data.get_att("lost_axis"),
      gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    gphys = gphys[true,-400..-1].rename("tr_v250").
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate. 
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)
  end

  def __gr_tr_om500
    # V
    gphys = @data.gphys_open("OMG").cut(true,0.5,true)
    lost_axis = [gphys.data.get_att("lost_axis"),
      gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    gphys = gphys[true,-400..-1].rename("tr_om500").
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate. 
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)
  end

  def __gr_spct(gphys)
    # 時空間スペクトル
    gphys = gphys.stspct_fft("#{gphys.name}_spct").
      set_att("ape_name",
	      "#{gphys.data.get_att("ape_name")}_S-T_power_spectrum").
      set_att("units","(#{gphys.data.get_att("units")} day-1)2").
      add_lost_axes(gphys.lost_axes)
    dim1 = gphys.coord(0).shape.to_s.to_i
    dim2 = gphys.coord(1).shape.to_s.to_i 
    mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..50])
  end

end