#!/usr/bin/env ruby

# ----------------------------------------------
# local load path

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

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

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

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


END{

  a = Ape_mkfig.new 3

    $rezol = "T39L48_eml"
    set_path
#    a.gr_comp_point("line")
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")
#    a.gr_comp_point
#    a.gr_comp_tuv

  a.gr_comp_tuv_horizontal("line")
  a.gr_comp_tuv_horizontal


  (1..7).each{ |num|

#    a.gr_comp_point("line#{num}")
#    a.gr_comp_tuv("line#{num}")
#    a.gr_comp_tuv("line#{num}_large")

    a.gr_comp_tuv_horizontal("line#{num}")

  }

  $rezol = "T39L48_ias"
  set_path
#  a.gr_comp_point("line")
#  a.gr_comp_tuv("line")
#  a.gr_comp_tuv("line_large") 
#  a.gr_comp_point
#  a.gr_comp_tuv
  
  a.gr_comp_tuv_horizontal("line")
  a.gr_comp_tuv_horizontal

  (1..3).each{ |num|
    
#    a.gr_comp_point("line#{num}")
#    a.gr_comp_tuv("line#{num}")
#    a.gr_comp_tuv("line#{num}_large")

    a.gr_comp_tuv_horizontal("line#{num}")

  }

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

    $rezol = "T39L48_#{item}"
    set_path
#    a.gr_comp_point
#    a.gr_comp_tuv
#    a.gr_comp_point("line")
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")

    a.gr_comp_tuv_horizontal("line")
    a.gr_comp_tuv_horizontal

  }



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

    $rezol = "T79L48_#{item}"
    set_path
#    a.gr_comp_point
#    a.gr_comp_tuv
#    a.gr_comp_point("line")
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")

    a.gr_comp_tuv_horizontal("line")
    a.gr_comp_tuv_horizontal


    $rezol = "T39L96_#{item}"
    set_path
#    a.gr_comp_point
#    a.gr_comp_tuv
#    a.gr_comp_point("line")
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")

    a.gr_comp_tuv_horizontal("line")
    a.gr_comp_tuv_horizontal

    $rezol = "T39L24_#{item}"
    set_path
#    a.gr_comp_point
#    a.gr_comp_tuv
#    a.gr_comp_point("line")
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")

    a.gr_comp_tuv_horizontal("line")
    a.gr_comp_tuv_horizontal

    $rezol = "T159L48_#{item}"
    set_path
#    a.gr_comp_point
#    a.gr_comp_point("line")
#    a.gr_comp_tuv
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")

  }


  $rezol = "T159L48_non"
  set_path
#    a.gr_comp_point
#    a.gr_comp_tuv
#    a.gr_comp_point("line")
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")

  a.gr_comp_tuv_horizontal("line")
  a.gr_comp_tuv_horizontal

    $rezol = "T159L48_eml"
    set_path
#    a.gr_comp_point("line")
#    a.gr_comp_tuv("line")
#    a.gr_comp_tuv("line_large")

  a.gr_comp_tuv_horizontal("line")
  a.gr_comp_tuv_horizontal


 ["T39L96_non","T159L48_eml","T159L48_non"].each{ |rezol|
#  ["T159L48_eml","T159L48_non"].each{ |rezol|
# ["T39L96_non"].each{ |rezol|
    
    $rezol = rezol
    set_path
#    a.gr_comp_point("line1")
#    a.gr_comp_tuv("line1")
#    a.gr_comp_tuv("line1_large")

  a.gr_comp_tuv_horizontal("line1")

  }

}

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


def set_path
  $expID = "control"
  $fig_path = "/home/yukiko/work/ape/yukiko/figs/tmp/"
  $gr2ncfile_path = "/home/yukiko/work/ape/yukiko/data/#{$rezol}_expID01/"
  $compfile_path = "/home/yukiko/work/ape/yukiko/data/composite/"
  $file_label = "#{$rezol}_#{$expID}"
  $groupid = "AGUforAPE-03a"
end

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

class Ape_mkfig

  def mkfig_plot_comp(gphys, gphys_comp)

    if @fig_num == 3 then
      @data.mkdump_comp(gphys, gphys_comp)
    else
      @data.plot_comp(gphys, gphys_comp)
    end

  end

  def gr_comp_point(compmethod=nil)
    
    if compmethod == nil 
      compmethod="threshold"
    end
   
    if $rezol == "T159L48_eml"
      rain = __gphys_opn_marge("tr_tppn",0)[true,-801..-401]
    else
      rain = __gphys_opn_marge("tr_tppn",0)[true,-401..-1]
    end
    lost_axes = rain.lost_axes.to_s.sub("y=","lat=")
    rain = rain.lon_lotate.
      rename("comp_point_tr_tppn_#{compmethod}")

    # tr_u250
    if $rezol == "T159L48_eml"
      u250 = __gphys_opn_marge("U",0.25)[true,-801..-401]
    else
      u250 = __gphys_opn_marge("U",0.25)[true,-401..-1]
    end
    lost_axis = [u250.data.get_att("lost_axis"),
      u250.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    u250 = ( u250 - u250.mean(0) ).lon_lotate.
      add_lost_axes("(diff) from (mean) zonal"). 
      rename("comp_point_tr_u250_#{compmethod}")

    # tr_u850
    if $rezol == "T159L48_eml"
      u850 = __gphys_opn_marge("U",0.85)[true,-801..-401]
    else
      u850 = __gphys_opn_marge("U",0.85)[true,-401..-1]
    end
    lost_axis = [u850.data.get_att("lost_axis"),
      u850.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    u850 = ( u850 - u850.mean(0) ).lon_lotate.
      add_lost_axes("(diff) from (mean) zonal"). 
      rename("comp_point_tr_u850_#{compmethod}")

    # tr_q850
    if $rezol == "T159L48_eml"
      q850 = __gphys_opn_marge("Q",0.85)[true,-801..-401]
    else
      q850 = __gphys_opn_marge("Q",0.85)[true,-401..-1]
    end
    lost_axis = [q850.data.get_att("lost_axis"),
      q850.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    q850 = ( q850 - q850.mean(0) ).lon_lotate.
      add_lost_axes("(diff) from (mean) zonal"). 
      rename("comp_point_tr_q850_#{compmethod}")

    # tr_t570
    if $rezol == "T159L48_eml"
      t570 = __gphys_opn_marge("T",0.57)[true,-801..-401]
    else
      t570 = __gphys_opn_marge("T",0.57)[true,-401..-1]
    end
    lost_axis = [t570.data.get_att("lost_axis"),
      t570.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    t570 = ( t570 - t570.mean(0) ).lon_lotate.
      add_lost_axes("(diff) from (mean) zonal"). 
      rename("comp_point_tr_t570_#{compmethod}")



#------------------------------------------------------------------
=begin
    # tr_u350
    u350 = __gphys_opn_marge("U").cut(true,0.35,true).lon_lotate
    lost_axis = [u350.data.get_att("lost_axis"),
      u350.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    u350 = u350[true,-401..-1].rename("tr_u350")
    u350 = ( u350 - u350.mean(0) )
    grid_1 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time").
            put_att("units","days"))
    grid_0 = u350.grid_copy.axis(0)
    grid = Grid.new(grid_0,grid_1)
    u350 = GPhys.new(grid, u350.data ).set_lost_axes(lost_axis).
      add_lost_axes("(diff) from (mean) zonal"). 
      rename("comp_point_tr_u350_#{compmethod}")


    # tr_t350
    t350 = __gphys_opn_marge("T").cut(true,0.35,true).lon_lotate
    lost_axis = [t350.data.get_att("lost_axis"),
      t350.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    t350 = t350[true,-401..-1].rename("tr_t350")
    t350 = ( t350 - t350.mean(0) )
    grid_1 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time").
            put_att("units","days"))
    grid_0 = t350.grid_copy.axis(0)
    grid = Grid.new(grid_0,grid_1)
    t350 = GPhys.new(grid, t350.data ).set_lost_axes(lost_axis).
      add_lost_axes("(diff) from (mean) zonal"). 
      rename("comp_point_tr_t350_#{compmethod}")
=end
#------------------------------------------------------------------


    file_name = 
      "#{$compfile_path}/#{$rezol}_expID01_composite_#{compmethod.gsub("_large","")}.nc"
    @data = Ape.new(file_name)
    comp = @data.gphys_open("comp_point")

    lost_axis = comp.data.get_att("lost_axis").split(/\n/)[1]
    
    rain = rain.add_lost_axes(lost_axis )
    u250 = u250.add_lost_axes(lost_axis )
    u850 = u850.add_lost_axes(lost_axis )
    q850 = q850.add_lost_axes(lost_axis )
    t570 = t570.add_lost_axes(lost_axis )

#    u350 = u350.add_lost_axes(lost_axis )
#    t350 = t350.add_lost_axes(lost_axis )


#    du850dx = u850 * 0.0 + 
#    ( u850.cshift(1).data.val - u850.cshift(-1).data.val )
#    du850dx = du850dx  / ( 2*( 40000*100 ) / grid_0.pos.val.size )
#    du850dx = du850dx.set_att("ape_name", "d(eastword_wind)/d(lon)"). 
#	  set_att("units","s-1").
#          rename("comp_point_tr_du850dx_#{compmethod}")

    mkfig_plot_comp(rain, comp)
    mkfig_plot_comp(u250, comp)
    mkfig_plot_comp(u850, comp)
    mkfig_plot_comp(q850, comp)
    mkfig_plot_comp(t570, comp)

#    mkfig_plot_comp(u350, comp)
#    mkfig_plot_comp(t350, comp)
#    mkfig_plot_comp(du850dx, comp)

  end


  def gr_comp_tuv(compmethod=nil)
    
    if compmethod == nil 
      compmethod="threshold"
    end

    file_name = 
      "#{$compfile_path}/#{$rezol}_expID01_composite_#{compmethod.gsub("_large","")}.nc"

    @data = Ape.new(file_name)
    t_comp = @data.gphys_open("comp_t").rename("comp_tuom_#{compmethod}")
    tv_comp = @data.gphys_open("comp_tv").rename("comp_tvuom_#{compmethod}")
    q_comp = @data.gphys_open("comp_q").rename("comp_quom_#{compmethod}")
    tconv_comp = @data.gphys_open("comp_tconv").
      rename("comp_tconvuom_#{compmethod}") 
    u_comp = @data.gphys_open("comp_u").rename("comp_uuom_#{compmethod}")
    om_comp = @data.gphys_open("comp_om")

    lost_axes = t_comp.get_att("lost_axis").split(/\n/)

    t_comp = t_comp.set_lost_axes(lost_axes).
      add_lost_axes("(diff) from (mean) zonal")
    tv_comp = tv_comp.set_lost_axes(lost_axes).
      add_lost_axes("(diff) from (mean) zonal")
    q_comp = q_comp.set_lost_axes(lost_axes).
      add_lost_axes("(diff) from (mean) zonal")
    u_comp = u_comp.set_lost_axes(lost_axes).
      add_lost_axes("(diff) from (mean) zonal")
    tconv_comp = tconv_comp.set_lost_axes(lost_axes)

    # ベクトルの間引
    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)
    
    if compmethod =~ /_large/ 
      num_l = u_comp.axis(0).pos.val.size/60
    else
      num_l = u_comp.axis(0).pos.val.size/120
    end
    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
    }

    # 描画
    dim = t_comp.axis(0).pos.val.size

    if compmethod =~ /_large/
      mkfig_plot(t_comp,
	         u_data,
	         -om_data )

      mkfig_plot(u_comp,
	         u_data,
	         -om_data )

      mkfig_plot(tv_comp,
	         u_data,
	         -om_data )

      mkfig_plot(q_comp, 
	         u_data, 
	         -om_data )

      mkfig_plot(tconv_comp, 
	         u_data, 
	         -om_data )
    else
      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])

      mkfig_plot(u_comp[(dim/3)..(dim*2/3),true],
	         u_data[(dim/3)..(dim*2/3),true], 
	         -om_data[(dim/3)..(dim*2/3),true])

      mkfig_plot(tv_comp[(dim/3)..(dim*2/3),true],
	         u_data[(dim/3)..(dim*2/3),true], 
	         -om_data[(dim/3)..(dim*2/3),true])

      mkfig_plot(q_comp[(dim/3)..(dim*2/3),true],
	         u_data[(dim/3)..(dim*2/3),true], 
	         -om_data[(dim/3)..(dim*2/3),true])

      mkfig_plot(tconv_comp[(dim/3)..(dim*2/3),true],
	         u_data[(dim/3)..(dim*2/3),true], 
	         -om_data[(dim/3)..(dim*2/3),true])
    end

  end

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

  def gr_comp_tuv_horizontal(compmethod=nil)
    
    if compmethod == nil 
      compmethod="threshold"
    end

    file_name = 
      "#{$compfile_path}/#{$rezol}_expID01_composite_horizontal_#{compmethod}.nc"

    @data = Ape.new(file_name)
    tppn_comp = @data.gphys_open("comp_tppn").
      rename("comp_tppn_#{compmethod}")
    psuv_comp = @data.gphys_open("comp_mslp").
      rename("comp_psuv_#{compmethod}")
    phiuv850_comp = @data.gphys_open("comp_phi850s").
      rename("comp_phiuv850_#{compmethod}")
    phiuv250_comp = @data.gphys_open("comp_phi250s").
      rename("comp_phiuv250_#{compmethod}")
    tuv570_comp = @data.gphys_open("comp_t570s").
      rename("comp_tuv570_#{compmethod}")
    quv850_comp = @data.gphys_open("comp_q850s").
      rename("comp_quv850_#{compmethod}")

    u250s_comp = @data.gphys_open("comp_u250s")
    u570s_comp = @data.gphys_open("comp_u570s")
    u850s_comp = @data.gphys_open("comp_u850s")
    u950s_comp = @data.gphys_open("comp_u950s")
    v250s_comp = @data.gphys_open("comp_v250s")
    v570s_comp = @data.gphys_open("comp_v570s")
    v850s_comp = @data.gphys_open("comp_v850s")
    v950s_comp = @data.gphys_open("comp_v950s")

    tppn_comp = tppn_comp.
      set_lost_axes( tppn_comp.get_att("lost_axis").split(/\n/) )
    psuv_comp = psuv_comp.
      set_lost_axes( psuv_comp.get_att("lost_axis").split(/\n/) )
    phiuv850_comp = phiuv850_comp.
      set_lost_axes( phiuv850_comp.get_att("lost_axis").split(/\n/) )
    phiuv250_comp = phiuv250_comp.
      set_lost_axes( phiuv250_comp.get_att("lost_axis").split(/\n/) )
    tuv570_comp = tuv570_comp.
      set_lost_axes( tuv570_comp.get_att("lost_axis").split(/\n/) )
    quv850_comp = quv850_comp.
      set_lost_axes( quv850_comp.get_att("lost_axis").split(/\n/) )

    # ベクトルの間引
    u250s_data = NArray.sfloat( u250s_comp.axis(0).pos.val.size , 
                                u250s_comp.axis(1).pos.val.size ).fill!(0.0)
    v250s_data = NArray.sfloat( v250s_comp.axis(0).pos.val.size , 
                                v250s_comp.axis(1).pos.val.size ).fill!(0.0)
    
    u850s_data = NArray.sfloat( u850s_comp.axis(0).pos.val.size , 
                                u850s_comp.axis(1).pos.val.size ).fill!(0.0)
    v850s_data = NArray.sfloat( v850s_comp.axis(0).pos.val.size , 
                                v850s_comp.axis(1).pos.val.size ).fill!(0.0)
    
    u950s_data = NArray.sfloat( u950s_comp.axis(0).pos.val.size , 
                                u950s_comp.axis(1).pos.val.size ).fill!(0.0)
    v950s_data = NArray.sfloat( v950s_comp.axis(0).pos.val.size , 
                                v950s_comp.axis(1).pos.val.size ).fill!(0.0)

    u570s_data = NArray.sfloat( u570s_comp.axis(0).pos.val.size , 
                                u570s_comp.axis(1).pos.val.size ).fill!(0.0)
    v570s_data = NArray.sfloat( v570s_comp.axis(0).pos.val.size , 
                                v570s_comp.axis(1).pos.val.size ).fill!(0.0)

    num_l = u250s_comp.axis(0).pos.val.size/120
    num_m = u250s_comp.axis(1).pos.val.size/60*3
    
    u250s_comp.axis(0).pos.val.size.times{ |num|
      if (num % num_l) == 0
	u250s_data[num,true]  =  u250s_comp.data.val[num,true]
	v250s_data[num,true]  =  v250s_comp.data.val[num,true]
	u850s_data[num,true]  =  u850s_comp.data.val[num,true]
	v850s_data[num,true]  =  v850s_comp.data.val[num,true]
	u950s_data[num,true]  =  u950s_comp.data.val[num,true]
	v950s_data[num,true]  =  v950s_comp.data.val[num,true]
	u570s_data[num,true]  =  u570s_comp.data.val[num,true]
	v570s_data[num,true]  =  v570s_comp.data.val[num,true]

      elsif
	u250s_data[num,true]  = 0.0
	v250s_data[num,true]  = 0.0
	u850s_data[num,true]  = 0.0
	v850s_data[num,true]  = 0.0
	u950s_data[num,true]  = 0.0
	v950s_data[num,true]  = 0.0
	u570s_data[num,true]  = 0.0
	v570s_data[num,true]  = 0.0
      end
    }
    
    num_m =1 if num_m < 1
    
    u250s_comp.axis(1).pos.val.size.times{ |num|
      if (num % num_m) == 0
	u250s_data[true,num]  =  u250s_data[true,num]
	v250s_data[true,num]  =  v250s_data[true,num]
	u850s_data[true,num]  =  u850s_data[true,num]
	v850s_data[true,num]  =  v850s_data[true,num]
	u950s_data[true,num]  =  u950s_data[true,num]
	v950s_data[true,num]  =  v950s_data[true,num]
	u570s_data[true,num]  =  u570s_data[true,num]
	v570s_data[true,num]  =  v570s_data[true,num]
      elsif
	u250s_data[true,num]  = 0.0
	v250s_data[true,num]  = 0.0
	u850s_data[true,num]  = 0.0
	v850s_data[true,num]  = 0.0
	u950s_data[true,num]  = 0.0
	v950s_data[true,num]  = 0.0
	u570s_data[true,num]  = 0.0
	v570s_data[true,num]  = 0.0
      end
    }

    # 描画
    dim = tppn_comp.axis(0).pos.val.size

    mkfig_plot(tppn_comp[(dim/3)..(dim*2/3),true] )

    mkfig_plot(psuv_comp[(dim/3)..(dim*2/3),true],
               u950s_data[(dim/3)..(dim*2/3),true], 
               v950s_data[(dim/3)..(dim*2/3),true])

    mkfig_plot(phiuv850_comp[(dim/3)..(dim*2/3),true],
               u850s_data[(dim/3)..(dim*2/3),true], 
               v850s_data[(dim/3)..(dim*2/3),true])

    mkfig_plot(phiuv250_comp[(dim/3)..(dim*2/3),true],
               u250s_data[(dim/3)..(dim*2/3),true], 
               v250s_data[(dim/3)..(dim*2/3),true])

    mkfig_plot(tuv570_comp[(dim/3)..(dim*2/3),true],
               u570s_data[(dim/3)..(dim*2/3),true], 
               v570s_data[(dim/3)..(dim*2/3),true])

    mkfig_plot(quv850_comp[(dim/3)..(dim*2/3),true],
               u850s_data[(dim/3)..(dim*2/3),true], 
               v850s_data[(dim/3)..(dim*2/3),true])

  end

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


  def __gphys_opn_marge(name,dimnum)

    if $rezol =~ /T159L48/

      g = 
        GPhys::NetCDF_IO.open(
                              "#{$gr2ncfile_path}#{$rezol}_expID01_0014.nc",
                             name).cut(true,dimnum,true)

      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)

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

      g = 
        GPhys::NetCDF_IO.open(
                              "#{$gr2ncfile_path}#{$rezol}_expID01_0007.nc",
                              name).cut(true,dimnum,true)

      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)

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

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

    return gphys

  end

end


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


class Ape
  
  def plot_comp(gphys, gphys_u=nil, gphys_v=nil)
    gropn(1) if $gropn == nil
    
    if gphys.class == Array
      plot_set(gphys[0]) 
    else 
      plot_set(gphys) 
    end

    plot_main_comp(gphys, gphys_u)

    # タイトルを標準出力へ
    if gphys.class == Array then
      if  gphys[0].data.get_att("ape_name") then
	puts "#{gphys[0].name} (#{gphys[0].data.get_att("ape_name")})"
      end
    elsif gphys.data.get_att("ape_name") then
      puts "#{gphys.name} (#{gphys.data.get_att("ape_name")})"
    else
      puts "#{gphys.name}"
    end
  end
  
  def mkdump_comp(gphys, gphys_u=nil, gphys_v=nil)
    gropn(3) if $gropn == nil

    plot_comp(gphys, gphys_u)

    gphys = gphys[0]     if gphys.class == Array
    
    file_name = "#{$file_label}_#{gphys.name}.gif"
    `rm dcl_* tmp.pnm`  
    plot(gphys)
    `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm;done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm`  

    $file_name = $pre_file

  end

  def plot_main_comp(gphys, gphys_u=nil)

    # attribute の long_name を消す (タイトル描画位置の都合)
    if gphys.data.get_att("long_name")
      gphys = gphys.set_att("long_name","")
    end

    GGraph.set_fig($fig_set_hash)

    # 描画
    GGraph.tone( gphys, true , $tone_hash ) 

    # タイトル
    DCL::uzinit 
    tani = "(#{gphys.data.get_att("units")})"
    DCL::uxsttl("t", tani , -1.0 )
    if  gphys.data.get_att("ape_name") then
      DCL::uxmttl("t", gphys.data.get_att("ape_name").gsub("_", " "), -1.0 ) 
    end

#    plot_set(gphys_u) 
    GGraph.contour( gphys_u, false , 
                 "lev"=>[0.5,0.8], 'label'=> ["",""], "index"=>5, 'nozero'=> true )
    
    # nc ファイル名
    $file_label = "#{@file}@#{gphys.name}"  if $file_label == "filename"
    DCL::uzinit
    DCL.uxpttl("t", 0, $file_label, 1.0)

    # トーンバー
    plot_set(gphys) 
    GGraph.color_bar($colorbar_hash)
#    DCL::Util::color_bar($cbar_conf)  

  end
end

class GPhys
  def cshift(n) 
    y = self.dup
    y[n..-1,true] = self[0..-1-n,true]
    y[0..n-1,true] = self[-n..-1,true]
    return y
  end
end

