#!/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
#    nc_tr_300day
#    gr_tr_all_300day

=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"  || item =~ "lon" || item =~ "lat" || item == "vertical" || item == "site" 

	gphys = @data.gphys_open(item)
	gphys = gphys[0,0,true] if $rezol == "K1JAPAN" || $rezol == "MGO"

	if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	  gphys = gphys[0,true] 
	end

	gphys = gphys[0,true] 	if $rezol == "UKMO_n96"


	if $rezol == "MGO" then 
	  grid_0 =
	    Axis.new().
	    set_pos(VArray.new(NArray.sfloat(1080).indgen!+180).rename("time").
		    put_att("units","days since 0000-01-01"))
	  grid = Grid.new(grid_0)
	  gphys = GPhys.new(grid, gphys.data )
	end

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	if $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "MGO" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96"
	  if item == "gt_sw_toar" || item == "gt_sswr" || item == "gt_slwd" 
	    mkfig_plot(-gphys)
	  else
	    mkfig_plot(gphys)
	  end
	else
	  mkfig_plot(gphys)
	end
      end
    }

    unless $rezol == "FRCGC" 
      # 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")
      gphys = gphys[0,0,true] if $rezol == "K1JAPAN" || $rezol == "MGO"
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[0,true] 
      end
      gphys = gphys[0,true] 	if $rezol == "UKMO_n96"

      if $rezol == "MGO" then 
	grid_0 =
	  Axis.new().
	  set_pos(VArray.new(NArray.sfloat(1080).indgen!+180).rename("time").
		  put_att("units","days since 0000-01-01"))
	grid = Grid.new(grid_0)
	gphys = GPhys.new(grid, gphys.data )
      end
      mkfig_plot(gphys)
    end
  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"  || item == "gw" || item =~ "time"  || item =~ "lev"  || item =~ "pressure"

	gphys = @data.gphys_open(item)
	gphys = gphys[true,true,true,0] if  $rezol == "MGO"

	if $rezol == "UKMO_n48" && $expID == "1keq" 
	  gphys = gphys[true,true,true,0] 
	end

	if $rezol == "UKMO_n96" 	    
	  gphys = gphys[true,true,true,0]  unless $expID == "control"
	end

	if ( $rezol == "CSIRO" || $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "NCAR" ) && item == "ml_rh"
	  gphys = gphys*100
	end
	if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96") && item == "ml_phi"
	  gphys = gphys*9.8
	end

	if $rezol == "FRCGC" && item == "ml_phi"
	  gphys = gphys*9.8
	end

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	mkfig_plot(gphys.mean(0))
      end
    }
    # mass_stream_function
    gphys = @data.gphys_open("ml_v")
    gphys = gphys[true,true,true,0] if $rezol == "MGO"
    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys = gphys[true,true,true,0] 
    end
    if $rezol == "UKMO_n96" 	    
      gphys = gphys[true,true,true,0]  unless $expID == "control"
    end
    gphys = gphys.mean(0).strm_function("ml_psi").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1") 
    if  $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ||  $rezol == "MGO"
      gphys = gphys*100
    end
    mkfig_plot(gphys)

    # MSE (moist static energy) = CpT + gz + Lq
    # DSE (dry static energy) = CpT + gz

    gphys_t = @data.gphys_open("ml_t").mean(0)
    gphys_q = @data.gphys_open("ml_q").mean(0)
    gphys_phi = @data.gphys_open("ml_phi").mean(0)

    if $rezol == "MGO"
      gphys_t = gphys_t[true,true,0] 
      gphys_q = gphys_q[true,true,0] 
      gphys_phi = gphys_phi[true,true,0] 
    end

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_t = gphys_t[true,true,0] 
      gphys_q = gphys_q[true,true,0] 
      gphys_phi = gphys_phi[true,true,0]   
    end
    
    if $rezol == "UKMO_n96" 	    
      unless $expID == "control"
	gphys_t = gphys_t[true,true,0] 
	gphys_q = gphys_q[true,true,0] 
	gphys_phi = gphys_phi[true,true,0]   
      end
    end

    if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96") 
      gphys_phi = gphys_phi*9.8
    end
    if $rezol == "FRCGC" && item == "ml_phi"
      gphys = gphys*9.8
    end

    # Cp: 1004 J kg-1 K-1
    # L : 2.500 * 10**6 J kg-1
    # g : 9.8 m s-2

    gphys_dse = gphys_phi + gphys_t.data.val*1004.0
    gphys_mse = gphys_dse + gphys_q.data.val*2.5*(10**6)

    mkfig_plot(gphys_dse.rename("ml_dse").
	       set_att("ape_name", "dry static energy").
	       set_att("units","J kg-1") )

    mkfig_plot(gphys_mse.rename("ml_mse").
	       set_att("ape_name", "moist static energy").
	       set_att("units","J kg-1") )
  end


  # #{$groupid}_TR_control.nc
  def nc_tr

    if $rezol == "UKMO_n96"
      @data = Ape.new(["#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc",
			"#{$ncfile_path}#{$groupid}_TR_#{$expID}_b.nc"])
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc")
    else
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")

      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
    end

#    @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")

    @data_tmp.netcdf_open.var_names.each { |name|
      if name =~ /lw_toa|mslp|500|250|ps/
	if name =~ /lw_toa|mslp|ps/
	  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)
	if $rezol == "FRCGC"
	  gphys =  gphys.cut(true,0,true)
	  grid_1 =
	    Axis.new().
	    set_pos(VArray.new(NArray.sfloat(120).indgen!/4).rename("time").
            put_att("units","days"))
	else
	  gphys =  gphys.cut(true,0,true)[true,-401..-1]
	  grid_1 =
	    Axis.new().
	    set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time").
            put_att("units","days"))
	end

	grid_0 = gphys.grid_copy.axis(0)
        grid = Grid.new(grid_0,grid_1)
        gphys = GPhys.new(grid, gphys.data )

	gphys = gphys.rename("tr_mslp") if gphys.name == "tr_ps" 

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	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
	if $rezol == "FRCGC"
	  mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..-36])
	else
	  mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
	end

      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
	if $rezol == "FRCGC"
	  gphys = gphys.cut(true,0,true)
	  grid_1 =
	    Axis.new().
	    set_pos(VArray.new(NArray.sfloat(120).indgen!/4).rename("time").
            put_att("units","days"))
	else
	  gphys = gphys.cut(true,0,true)[true,-401..-1]
     	  grid_1 =
	    Axis.new().
	    set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time").
            put_att("units","days"))
	end

	grid_0 = gphys.grid_copy.axis(0)	
        grid = Grid.new(grid_0,grid_1)
        gphys = GPhys.new(grid, gphys.data ) 
	
	if $rezol == "LASG" || $rezol == "NCAR"
	  gphys = gphys * 1000
	end

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	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
	if $rezol == "FRCGC"
	  puts dim2
	  mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..-36])
	else
	  mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
	end
      end
    }

  end


  # #{$groupid}_TR_control.nc
  def nc_tr_5deg_mean

    if $rezol == "UKMO_n96"
      @data = Ape.new(["#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc",
			"#{$ncfile_path}#{$groupid}_TR_#{$expID}_b.nc"])
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc")
    else
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
    end

    @data_tmp.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,-5..5,true)[true,true,-401..-1].
	  rename("#{name}_5deg_mean")
	lon = gphys.grid_copy.axis(0).pos.val 
	lat = gphys.grid_copy.axis(1).pos.val 
	lon.size.times { |lon_num| 
	    gphys[lon_num,true,true] = 
	      ( gphys[lon_num,true,true] * cos(lat*PI*2.0/360.0) )
	}
	gphys = gphys.mean(1)
	
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	if gphys.name == "tr_ps_5deg_mean" 
	  gphys = gphys.rename("tr_mslp_5deg_mean") 
	end

	mkfig_plot((gphys - gphys.mean(0)).
		   add_lost_axes(lost_axis).
		   add_lost_axes("(diff) from (mean) zonal").lon_lotate)

	# 時空間スペクトル
	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).
	  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..80])

      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
	gphys =  gphys.cut(true,-5..5,true)[true,true,-401..-1].
	  rename("#{name}_5deg_mean")
	lon = gphys.grid_copy.axis(0).pos.val 
	lat = gphys.grid_copy.axis(1).pos.val 
	lon.size.times { |lon_num| 
	  gphys[lon_num,true,true] = 
	    ( gphys[lon_num,true,true] * cos(lat*PI*2.0/360.0) )
	}
	gphys = gphys.mean(1)

	if $rezol == "LASG" || $rezol == "NCAR"
	  gphys = gphys * 1000
	end

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	mkfig_plot(gphys.lon_lotate)

	# 時空間スペクトル
	gphys = gphys.stspct_fft("tr_tppn_5deg_mean_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..80])
      end
    }

  end

  # #{$groupid}_TR_control.nc
  def nc_tr_300day

    if $rezol == "UKMO_n96"
      @data = Ape.new(["#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc",
			"#{$ncfile_path}#{$groupid}_TR_#{$expID}_b.nc"])
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc")
    else
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")

      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
    end

#    if $rezol == "T159L48_non"
#      @data = Ape.new(["#{$ncfile_path}#{$groupid}_TR_#{$expID}_half1.nc",
#			"#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc"])
#      data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half1.nc")
#    else 
#      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
#      data = @data
#    end

      @data_tmp.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,-801..-401].
	  rename("#{name}_200")
	gphys = gphys.rename("tr_mslp_200")  if gphys.name == "tr_ps_200" 
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	mkfig_plot((gphys - gphys.mean(0)).
		   add_lost_axes(lost_axis).
		   add_lost_axes("(diff) from (mean) zonal").lon_lotate)
	# 時空間スペクトル
	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).
	  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..80])

	#
	gphys = @data.gphys_open(name)
	gphys =  gphys.cut(true,0,true)[true,-1201..-801].
	  rename("#{name}_300")
	gphys = gphys.rename("tr_mslp_300")  if gphys.name == "tr_ps_300" 
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	mkfig_plot((gphys - gphys.mean(0)).
		   add_lost_axes(lost_axis).
		   add_lost_axes("(diff) from (mean) zonal").lon_lotate)
	# 時空間スペクトル
	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).
	  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..80])

      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
	gphys = gphys.cut(true,0,true)[true,-801..-401].rename("#{name}_200")
	if $rezol == "LASG" || $rezol == "NCAR"
	  gphys = gphys * 1000
	end
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	mkfig_plot(gphys.lon_lotate)
	# 時空間スペクトル
	gphys = gphys.stspct_fft("tr_tppn_200_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..80])

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
	gphys = gphys.cut(true,0,true)[true,-1201..-801].rename("#{name}_300")
	if $rezol == "LASG" || $rezol == "NCAR"
	  gphys = gphys * 1000
	end
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	mkfig_plot(gphys.lon_lotate)
	# 時空間スペクトル
	gphys = gphys.stspct_fft("tr_tppn_300_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..80])

      end
    }

  end


  # #{$groupid}_TR_control.nc
  def nc_tr_360day
    
    if $rezol == "UKMO_n96"
      @data = Ape.new(["#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc",
			"#{$ncfile_path}#{$groupid}_TR_#{$expID}_b.nc"])
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_a.nc")
    else
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
      
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
    end

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

#    data.netcdf_open.var_names.each { |name|
    @data_tmp.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).
	  rename("#{name}_360all")

	gphys = gphys.rename("tr_mslp_360all")  if gphys.name == "tr_ps_360all" 
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	mkfig_plot((gphys - gphys.mean(0)).
		   add_lost_axes(lost_axis).
		   add_lost_axes("(diff) from (mean) zonal").lon_lotate)

	# 時空間スペクトル
	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).
	  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..288])

      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
	gphys = gphys.cut(true,0,true).
	  rename("#{name}_360all")
	if $rezol == "LASG" || $rezol == "NCAR"
	  gphys = gphys * 1000
	end
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	mkfig_plot(gphys.lon_lotate)
	# 時空間スペクトル
	gphys = gphys.stspct_fft("tr_tppn_360all_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..288])
      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"  || item =~ "time"  || item == "nv"
	gphys = @data.gphys_open(item)
	gphys = gphys[true,true,0] if $rezol == "MGO"

	if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	  gphys = gphys[true,true,0] 
	end
	gphys = gphys[true,true,0] 	if $rezol == "UKMO_n96"

	gphys = gphys.lon_lotate
	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "MGO" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96"
	  if item == "sh_sw_toar" || item == "sh_sswr" || item == "sh_slwd" 
	      mkfig_plot(-gphys)
	    else
	      mkfig_plot(gphys)
	    end
	else
	      mkfig_plot(gphys)
	end

      end
    }
    unless $rezol == "FRCGC"
      # 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[true,true,0] if $rezol == "MGO"
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[true,true,0] 
      end
      gphys = gphys[true,true,0] 	if $rezol == "UKMO_n96"
      gphys = gphys.lon_lotate
      mkfig_plot(gphys)
    end
  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"  || item =~ "time"  || item == "gw"  || item == "nv" 
	puts item
	gphys = @data.gphys_open(item)
	if  $rezol == "MGO"
	  gphys = gphys[true,true,0] 
	end
	if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	  gphys = gphys[true,true,0] 
	end
	gphys = gphys[true,true,0] 	if $rezol == "UKMO_n96"

	gphys = gphys.rename("#{gphys.name}_zonal")

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	if $rezol == "CSIRO" || $rezol == "GSFC" || $rezol == "LASG" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96"
	  if item == "sh_tauv" || item == "sh_tauu"
	    gphys = -gphys
	  end
	end

	if $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "MGO" || $rezol == "NCAR" || $rezol == "UKMO_n96" || $rezol == "UKMO_n48"
	  if item == "sh_sw_toar" || item == "sh_sswr" || item == "sh_slwd" 
	      mkfig_plot(-gphys.mean(0))
	    else
	      mkfig_plot(gphys.mean(0))
	    end
	else
	      mkfig_plot(gphys.mean(0))
	end

      end
    }
    unless $rezol == "FRCGC"
      # 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")
      gphys = gphys[true,true,0] if $rezol == "MGO"
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[true,true,0] 
      end
      gphys = gphys[true,true,0] 	if $rezol == "UKMO_n96"
      mkfig_plot(gphys.mean(0))
    end
  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"  || item =~ "lon" || item =~ "lat"
# gt_slwu は一定値

	gphys = @data.gphys_open(item)
	gphys_con = @data_con.gphys_open(item)

	if $rezol == "K1KAPAN" || $rezol == "MGO"
	  gphys = gphys[0,0,true] ; gphys_con = gphys_con[0,0,true]
	end

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	gphys_con = gphys_con.mean(0)

	lost_axis = ["(diff) from (mean) time of #{zonalcont}"]
	gphys = gphys.add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys - gphys_con)

	if item == "gt_sw_toar" || item == "gt_sswr" || item == "gt_slwd" 
	  mkfig_plot(-gphys + gphys_con)
	end

      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("gt_cppn")
    g2 = @data.gphys_open("gt_dppn")

    if $rezol == "K1KAPAN" || $rezol == "MGO"
      g1 = g1[0,0,true] ; g2 = g2[0,0,true]
    end

    gphys = (g1 + g2)
    g1 = @data_con.gphys_open("gt_cppn")
    g2 = @data_con.gphys_open("gt_dppn")

    if $rezol == "K1KAPAN" || $rezol == "MGO"
      g1 = g1[0,0,true] ; g2 = g2[0,0,true]
    end

    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"  || item == "gw" || item =~ "time"  || item =~ "lev"   || item =~ "pressure"
	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")

	if $rezol == "UKMO_n48" && $expID == "1keq" 
	  gphys = gphys[true,true,true,0] 
	end
    
	if $rezol == "UKMO_n96" 	    
	  gphys = gphys[true,true,true,0] 	  unless $expID == "control"
	end

	if ( $rezol == "CSIRO" || $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "NCAR" ) && item == "ml_rh"
	  gphys = gphys*100
	end
	if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ) && item == "ml_phi"
	  gphys = gphys*9.8
	  gphys_con = gphys_con*9.8
	end

	if $rezol == "FRCGC" && item == "ml_phi"
	  gphys = gphys*9.8
	  gphys_con = gphys_con*9.8
	end

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	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")
    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys = gphys[true,true,true,0] 
    end
    if $rezol == "UKMO_n96" 	    
      gphys = gphys[true,true,true,0] 	  unless $expID == "control"
    end

    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)

    # MSE (moist static energy) = CpT + gz + Lq
    # DSE (dry static energy) = CpT + gz

    gphys_t = @data.gphys_open("ml_t").mean(0)
    gphys_q = @data.gphys_open("ml_q").mean(0)
    gphys_phi = @data.gphys_open("ml_phi").mean(0)

    gphys_t_con = @data_con.gphys_open("ml_t").mean(0)
    gphys_q_con = @data_con.gphys_open("ml_q").mean(0)
    gphys_phi_con = @data_con.gphys_open("ml_phi").mean(0)

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_t = gphys_t[true,true,0] 
      gphys_q = gphys_q[true,true,0] 
      gphys_phi = gphys_phi[true,true,0]   
    end
    
    if $rezol == "UKMO_n96" 	    
      unless $expID == "control"
	gphys_t = gphys_t[true,true,0] 
	gphys_q = gphys_q[true,true,0] 
	gphys_phi = gphys_phi[true,true,0]   
      end
    end

    if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96") 
      gphys_phi = gphys_phi*9.8
      gphys_phi_con = gphys_phi_con*9.8
    end
    if $rezol == "FRCGC" && item == "ml_phi"
      gphys = gphys*9.8
      gphys_phi_con = gphys_phi_con*9.8
    end


    # Cp: 1004 J kg-1 K-1
    # L : 2.500 * 10**6 J kg-1
    # g : 9.8 m s-2

    gphys_dse = gphys_phi + gphys_t.data.val*1004.0
    gphys_mse = gphys_dse + gphys_q.data.val*2.5*(10**6)

    gphys_dse_con = gphys_phi_con + gphys_t_con.data.val*1004.0
    gphys_mse_con = gphys_dse_con + gphys_q_con.data.val*2.5*(10**6)

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

    mkfig_plot(gphys_dse.rename("ml_dse_anm").
	       add_lost_axes(lost_axis).
	       set_att("ape_name", "dry static energy").
	       set_att("units","J kg-1") - gphys_dse_con)

    mkfig_plot(gphys_mse.rename("ml_mse_anm").
	       add_lost_axes(lost_axis).
	       set_att("ape_name", "moist static energy").
	       set_att("units","J kg-1") - gphys_mse_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"  || item =~ "time" || item == "nv"  
	gphys_con = @data_con.gphys_open(item)
	gphys_con = gphys_con[true,true,0] if $rezol == "MGO"
	lost_axis = ["(diff) from (mean) zonal of #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys[true,true,0] if $rezol == "MGO"
	if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	  gphys = gphys[true,true,0] 
	end
	if $rezol == "UKMO_n96"
	  gphys = gphys[true,true,0]
	  gphys_con = gphys_con[true,true,0]
	end

	gphys = (gphys - gphys_con.mean(0))
	gphys = gphys.lon_lotate.rename("#{gphys.name}_anm")

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	if $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "MGO" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96"
	  if item == "sh_sw_toar" || item == "sh_sswr" || item == "sh_slwd" 
	      mkfig_plot(-gphys)
	    else
	      mkfig_plot(gphys)
	    end
	else
	      mkfig_plot(gphys)
	end

      end
    }
    unless $rezol == "FRCGC"
      # total_precipitation_flux
      g1 = @data.gphys_open("sh_cppn")
      g2 = @data.gphys_open("sh_dppn")
      gphys = (g1 + g2)
      gphys = gphys[true,true,0] if $rezol == "MGO"
      g1 = @data_con.gphys_open("sh_cppn")
      g2 = @data_con.gphys_open("sh_dppn")
      gphys_con = (g1 + g2)
      gphys_con = gphys_con[true,true,0] if $rezol == "MGO"
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[true,true,0] 
      end
      if $rezol == "UKMO_n96"
	gphys = gphys[true,true,0]
	gphys_con = gphys_con[true,true,0]
      end

     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 = gphys[true,true,0] if $rezol == "MGO"
      gphys_con = @data_con.gphys_open("sh_slh")
      gphys_con = gphys_con[true,true,0] if $rezol == "MGO"
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[true,true,0] 
      end
      if $rezol == "UKMO_n96"
	gphys = gphys[true,true,0]
	gphys_con = gphys_con[true,true,0]
      end

      gphys = (gphys - gphys_con.mean(0)).lon_lotate
      if $rezol == "MGO"
	gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
      end
      if $rezol =~ /UKMO/
	gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
      end

      gphys = gphys.cut(true,0).rename("sh_slh_lat0_anm").
	add_lost_axes(lost_axis)
      mkfig_plot(gphys)
    end
  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"  || item =~ "time" || item =~  "bnds"  || item == "gw"  || item == "nv" #|| item =~ "sh_sw_toai"
	
	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")

	gphys = gphys[true,true,0] if $rezol == "MGO"
	gphys_con = gphys_con[true,true,0] if $rezol == "MGO"
	if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	  gphys = gphys[true,true,0] 
	end
	if $rezol == "UKMO_n96"
	  gphys = gphys[true,true,0]
	  gphys_con = gphys_con[true,true,0]
	end

	if $rezol == "MGO"
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end
	if $rezol =~ /UKMO/
	  gphys = gphys.set_att("ape_name", "#{gphys.data.get_att("long_name")}")
	end

	if $rezol == "CSIRO" || $rezol == "GSFC" || $rezol == "LASG" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96"
	  if item == "sh_tauv" || item == "sh_tauu"
	    gphys = -gphys
	    gphys_con = -gphys_con
	  end
	end

	if $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "MGO" || $rezol == "NCAR" || $rezol == "UKMO_n96" || $rezol == "UKMO_n48"
	  if item == "sh_sw_toar" || item == "sh_sswr" || item == "sh_slwd" 
	    mkfig_plot(- gphys.mean(0)+ gphys_con.mean(0))
	  else
	    mkfig_plot(gphys.mean(0)- gphys_con.mean(0))
	  end
	elsif $rezol =~ "non"
	  if item == "sh_cppn"
	  else 
	    mkfig_plot(gphys.mean(0)- gphys_con.mean(0)) 
	  end	    
	else
	  mkfig_plot(gphys.mean(0)- gphys_con.mean(0)) 
	end

      end
    }
    unless $rezol == "FRCGC"    

      # total_precipitation_flux
      if $rezol =~ "non"
	gphys = @data.gphys_open("sh_dppn")
	gphys_con = @data_con.gphys_open("sh_dppn")
      else 
	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)
      end

      gphys = gphys[true,true,0] if $rezol == "MGO"
      gphys_con = gphys_con[true,true,0] if $rezol == "MGO"
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[true,true,0] 
      end
      if $rezol == "UKMO_n96"
	gphys = gphys[true,true,0]
	gphys_con = gphys_con[true,true,0]
      end

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

    if $rezol == "MGO"
      gphys_ps = gphys_ps[true,true,0]
      gphys_ps_con = gphys_ps_con[true,true,0] 
    end
    if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
      gphys_ps = gphys_ps[true,true,0] 
    end

    if $rezol == "UKMO_n96" 	    
	gphys_ps = gphys_ps[true,true,0] 
	gphys_ps_con = gphys_ps_con[true,true,0] 
    end

    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, m s-1)")

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

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,100000).lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,100000).lon_lotate
    elsif $rezol == "UKMO_n96" 
      unless $expID == "control"
	gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	  cut(true,true,100000).lon_lotate
	gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	  cut(true,true,100000).lon_lotate
      end
    else
      gphys_u = @data.gphys_open("ml_u").cut(true,true,100000).lon_lotate
      gphys_v = @data.gphys_open("ml_v").cut(true,true,100000).lon_lotate
    end

    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,100000).lon_lotate
    gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,100000).lon_lotate


    # これでいいのかなぁ, 欠損値処理
    mask = gphys_u.gt(-900) * gphys_u.lt(10000)
    gphys_u = NArrayMiss.to_nam(gphys_u.data.val, mask)
    mask = gphys_u_con.gt(-900) * gphys_u_con.lt(10000)
    gphys_u_con = NArrayMiss.to_nam(gphys_u_con.data.val, mask)

    mask = gphys_v.gt(-900) * gphys_v.lt(10000)
    gphys_v = NArrayMiss.to_nam(gphys_v.data.val, mask)
    mask = gphys_v_con.gt(-900) * gphys_v_con.lt(10000)
    gphys_v_con = NArrayMiss.to_nam(gphys_v_con.data.val, mask)

    gphys_u_con[0,true].size.times{ |num|
      gphys_u_con[num,true] =  gphys_u_con.mean(0)
    }
    gphys_v_con[0,true].size.times{ |num|
      gphys_v_con[num,true] =  gphys_v_con.mean(0)
    }

    gphys_u = ( gphys_u - gphys_u_con )
    gphys_v = ( gphys_v - gphys_v_con )

    u_data  = NArrayMiss.sfloat( gphys_u[true,0].size ,
			    gphys_u[0,true].size ).fill!(0.0)
    v_data  = NArrayMiss.sfloat( gphys_v[true,0].size ,
			    gphys_v[0,true].size ).fill!(0.0)
    
    num_x = gphys_u[true,0].size/40
    num_y = gphys_u[0,true].size/20

    gphys_u[true,0].size.times{ |num|
      if (num % num_x) == 0
	u_data[num,true]  =  gphys_u[num,true]
	v_data[num,true]  =  gphys_v[num,true]
      elsif   
	u_data[num,true]  = 0.0
	v_data[num,true]  = 0.0
      end
    }  
    
    gphys_u[0,true].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
    }

    DCL::gllset("lmiss",true)
    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)

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].cut(true,true,25000)
    elsif $rezol == "UKMO_n96" 	    
      unless $expID == "control"
	gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].cut(true,true,25000)
      end
    else
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,25000)
    end

    if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ) 
      gphys_phi = gphys_phi*9.8
      gphys_phi_con = gphys_phi_con*9.8
    end

    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, m s-1)")

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,25000).lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,25000).lon_lotate
    elsif $rezol == "UKMO_n96" 	    
      unless $expID == "control"
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,25000).lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,25000).lon_lotate
      end
    else
      gphys_u = @data.gphys_open("ml_u").cut(true,true,25000).lon_lotate
      gphys_v = @data.gphys_open("ml_v").cut(true,true,25000).lon_lotate
    end

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

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


    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)

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].
	cut(true,true,85000)
    elsif $rezol == "UKMO_n96" 	    
      unless $expID == "control"
      gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].
	cut(true,true,85000)
      end
    else
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,85000)
    end

    if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ) 
      gphys_phi = gphys_phi*9.8
      gphys_phi_con = gphys_phi_con*9.8
    end

    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, m s-1)")

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,85000)
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,85000)
    elsif $rezol == "UKMO_n96" 	    
      unless $expID == "control"
	gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	  cut(true,true,85000)
	gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	  cut(true,true,85000)
      end
    else
      gphys_u = @data.gphys_open("ml_u").cut(true,true,85000)
      gphys_v = @data.gphys_open("ml_v").cut(true,true,85000)
    end

    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)

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_t = @data.gphys_open("ml_t")[true,true,true,0].
	cut(true,true,50000)
    elsif $rezol == "UKMO_n96" 	    
      unless $expID == "control"
      gphys_t = @data.gphys_open("ml_t")[true,true,true,0].
	cut(true,true,50000)
      end
    else
      gphys_t = @data.gphys_open("ml_t").cut(true,true,50000)
    end

    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, m s-1)")

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,50000)
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,50000)
    elsif $rezol == "UKMO_n96" 	    
      unless $expID == "control"
	gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	  cut(true,true,50000)
	gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	  cut(true,true,50000)
      end
    else
      gphys_u = @data.gphys_open("ml_u").cut(true,true,50000)
      gphys_v = @data.gphys_open("ml_v").cut(true,true,50000)
    end

    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")
    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_t = gphys_t[true,true,true,0] 
    end

    if $rezol == "UKMO_n96" 	    
      gphys_t = gphys_t[true,true,true,0]       unless $expID == "control"
    end
    
    gphys_t = gphys_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, Pa s-1)")

    gphys_u = @data.gphys_open("ml_u")
    gphys_om = @data.gphys_open("ml_om")
    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_u = gphys_u[true,true,true,0] 
      gphys_om = gphys_om[true,true,true,0] 
    end

    if $rezol == "UKMO_n96" 	    
      unless $expID == "control"
	gphys_u = gphys_u[true,true,true,0] 
	gphys_om = gphys_om[true,true,true,0] 
      end
    end

    gphys_u = gphys_u.cut(true,0,true)
    gphys_om = gphys_om.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_om").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","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1"]

    id_ary = Array.new

    id_ary.push( ["control","peaked","flat","Qobs","control-5N"])
    id_ary.push( ["control","1keq","3keq","3kw1"])
    if $rezol =~ /T39L48/
      id_ary.push( ["control","flat","Qobs","H1998con"])
      id_ary.push( ["control","3keq","flat","flat3keq","Qobs","Qobs3keq"]) 
      id_ary.push( ["control","3keq","H1998con","H1998pa"])
    end

    $file_label = "#{$rezol}"
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_control.nc")

    id_ary.size.times{ |idsize|
      @data.netcdf_open.var_names.each { |item| 
	id = id_ary[idsize]
	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).
			   rename("#{item}_#{idsize}").
			   set_att("line_name",$expID) )
	  }
	  puts item
	  if $rezol =~ /non/ && item == "gt_cppn"
	  else
	    mkfig_plot(gphys_ary) 
	  end
	end
      }
    }

    # total_precipitation_flux
    id_ary.size.times{ |idsize|
      id = id_ary[idsize]
      gphys_ary = []
      id.each{ |expID|
	$expID = expID
	@data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
	if $rezol =~ /non/
	  gphys = @data.gphys_open("gt_dppn").rename("gt_tppn")
	else
	  g1 = @data.gphys_open("gt_cppn")
	  g2 = @data.gphys_open("gt_dppn")
	  gphys = (g1 + g2).rename("gt_tppn")
	end
	gphys = gphys.
	  set_att("ape_name", "total_precipitation_flux")
	gphys_ary.push(gphys.
		       rename("#{gphys.name}_#{idsize}").
		       set_att("line_name",$expID) )
      }
      mkfig_plot(gphys_ary)
    } 
  end
  

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal_comparing_expid
#    id = ["control","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1"]
#    id = ["control","flat","peaked","Qobs","control-5N"]
#    id = ["control","flat","3keq","Qobs"]

    id_ary = Array.new

    id_ary.push( ["control","peaked","flat","Qobs","control-5N"])
    id_ary.push( ["control","1keq","3keq","3kw1"])
    if $rezol =~ /T39L48/
      id_ary.push( ["control","flat","Qobs","H1998con"])
      id_ary.push( ["control","3keq","flat","flat3keq","Qobs","Qobs3keq"]) 
      id_ary.push( ["control","3keq","H1998con","H1998pa"])
    end


    $file_label = "#{$rezol}"
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_control.nc")
    id_ary.size.times{ |idsize|
      @data.netcdf_open.var_names.each { |item| 
	id = id_ary[idsize]
	gphys_ary = []
	unless item =~ "lon" || item =~ "lat"  || item == "nv"
	  id.each{ |expID|
	    $expID = expID
	    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
	    gphys = @data.gphys_open(item).mean(0)
	    if  $rezol == "MGO"
	      gphys = gphys[true,0] 
	    end
	    if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	      gphys = gphys[true,0] 
	    end
	    gphys_ary.push(gphys.rename("#{gphys.name}_zonal").
			   rename("#{gphys.name}_#{idsize}").
			   set_att("line_name",$expID) )
	  }
	  puts item
	  if $rezol =~ /non/ && item == "sh_cppn"
	  else
	    mkfig_plot(gphys_ary)
	  end
	end

      }
    }
    # total_precipitation_flux
    id_ary.size.times{ |idsize|
      gphys_ary = []
      id = id_ary[idsize]
      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")
	if  $rezol == "MGO"
	  gphys = gphys[true,true,0] 
	end
	if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	  gphys = gphys[true,true,0] 
	end
	gphys_ary.push(gphys.mean(0).
		       rename("#{gphys.name}_#{idsize}").
		       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)

      if  $rezol == "MGO"
	gphys = gphys[true,true,0] 
	gphys_con = gphys_con[true,true,0] 
      end
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[true,true,0] 
      end

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

      if  $rezol == "MGO"
	gphys = gphys[true,true,0] 
	gphys_con = gphys_con[true,true,0] 
      end
      if $rezol == "UKMO_n48" && ($expID == "flat" ||$expID == "qobs" ||$expID == "1keq" ||$expID == "peaked")
	gphys = gphys[true,true,0] 
	gphys_con = gphys_con[true,true,0] 
      end

      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_ml_zonal_comparing_expid

    id = ["control","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1"]
    $file_label = "#{$rezol}"
    gphys_ary = [] 
    # U の赤道鉛直断面
    id.each{ |expID|
      $expID = expID
      $groupid = $groupid_hash[$rezol]
      $ncfile_path = "/work11/ape/NetCDF/#{$rezol}/"
      puts "#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc"
      @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",$expID) )
    }
    mkfig_plot(gphys_ary)

    # MSE (moist static energy) = CpT + gz + Lq
    # DSE (dry static energy) = CpT + gz
    gphys_dse_ary = []
    gphys_mse_ary = []
    id.each{ |expID|
      $expID = expID
      $groupid = $groupid_hash[$rezol]
      $ncfile_path = "/work11/ape/NetCDF/#{$rezol}/"
      puts "#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      gphys_t = @data.gphys_open("ml_t").mean(0).cut(20,true)
      gphys_q = @data.gphys_open("ml_q").mean(0).cut(20,true)
      gphys_phi = @data.gphys_open("ml_phi").mean(0).cut(20,true)
      if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ) 
	gphys_phi = gphys_phi*9.8
      end
      gphys_dse = gphys_phi + gphys_t.data.val*1004.0
      gphys_mse = gphys_dse + gphys_q.data.val*2.5*(10**6)
      gphys_dse_ary.push( gphys_dse.rename("ml_dse_lat20_zonal").
			 set_att("line_name",$expID).
			 set_att("ape_name", "dry static energy").
			 set_att("units","J kg-1") )
      gphys_mse_ary.push( gphys_mse.rename("ml_mse_lat20_zonal").
			 set_att("line_name",$expID).
			 set_att("ape_name", "moist static energy").
			 set_att("units","J kg-1") )
      puts $rezol
    }
    mkfig_plot(gphys_dse_ary)
    mkfig_plot(gphys_mse_ary)

  end


  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal_comparing_rezol
    $file_label = "#{$cumulus}_#{$expID}"
#    rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}"]
    rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}","T319L48_#{$cumulus}"]


    $ncfile_path = "/work11/ape/NetCDF/#{rezol[0]}/"

    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item =~ "lon" || item =~ "lat" || item == "nv"
	rezol.each{ |rezol_item|
	  $rezol = rezol_item
	  $ncfile_path = "/work11/ape/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 = "/work11/ape/NetCDF/#{$rezol}/"
#      $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_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}"]
    rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}","T319L48_#{$cumulus}"]

    $ncfile_path = "/work11/ape/NetCDF/#{rezol[0]}/"
    print "#{$ncfile_path}#{$groupid}\n"
    @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 = "/work11/ape/NetCDF/#{$rezol}/"
#	  $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 = "/work11/ape/NetCDF/#{$rezol}/"
#      $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_#{$cumulus}", "T39L24_#{$cumulus}","T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}"]


    # U の赤道鉛直断面
    gphys_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
#      $ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
      $ncfile_path = "/work11/ape/NetCDF/#{$rezol}/"
      puts $ncfile_path
      @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)

    # MSE (moist static energy) = CpT + gz + Lq
    # DSE (dry static energy) = CpT + gz
    gphys_dse_ary = []
    gphys_mse_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/work11/ape/NetCDF/#{$rezol}/"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      gphys_t = @data.gphys_open("ml_t").mean(0).cut(20,true)
      gphys_q = @data.gphys_open("ml_q").mean(0).cut(20,true)
      gphys_phi = @data.gphys_open("ml_phi").mean(0).cut(20,true)
      if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ) 
	gphys_phi = gphys_phi*9.8
      end
      gphys_dse = gphys_phi + gphys_t.data.val*1004.0
      gphys_mse = gphys_dse + gphys_q.data.val*2.5*(10**6)
      gphys_dse_ary.push( gphys_dse.rename("ml_dse_lat20_zonal").
			 set_att("line_name",$rezol).
			 set_att("ape_name", "dry static energy").
			 set_att("units","J kg-1") )
      gphys_mse_ary.push( gphys_mse.rename("ml_mse_lat20_zonal").
			 set_att("line_name",$rezol).
			 set_att("ape_name", "moist static energy").
			 set_att("units","J kg-1") )
      puts $rezol
    }
    mkfig_plot(gphys_dse_ary)
    mkfig_plot(gphys_mse_ary)

  end


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

    rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}","T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}"]

    expID_con_hash =  {
      "1keq" => "control", 
      "3keq" => "control", 
      "3kw1" => "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 = "/work11/ape/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_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" || item == "nv" 
	rezol.each{ |rezol_item|
	  $rezol = rezol_item
	  $ncfile_path = "/work11/ape/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 = "/work11/ape/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 = "/work11/ape/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 = "/work11/ape/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 =  {
      "1keq" => "control", 
      "3keq" => "control", 
      "3kw1" => "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 = "/work11/ape/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_ml_zonal_comparing_cumulus


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

    $file_label = "T39L48_#{$expID}"

    gphys_ary = [] 
    # U の赤道鉛直断面
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/work11/ape/NetCDF/#{$rezol}/"
      $groupid = $groupid_hash[$rezol]
      puts "#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc"
      @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)

    # MSE (moist static energy) = CpT + gz + Lq
    # DSE (dry static energy) = CpT + gz
    gphys_dse_ary = []
    gphys_mse_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/work11/ape/NetCDF/#{$rezol}/"
      $groupid = $groupid_hash[$rezol]
      puts "#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc"
      @data     = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      gphys_t = @data.gphys_open("ml_t").mean(0).cut(20,true)
      gphys_q = @data.gphys_open("ml_q").mean(0).cut(20,true)
      gphys_phi = @data.gphys_open("ml_phi").mean(0).cut(20,true)
      if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ) 
	gphys_phi = gphys_phi*9.8
      end
      gphys_dse = gphys_phi + gphys_t.data.val*1004.0
      gphys_mse = gphys_dse + gphys_q.data.val*2.5*(10**6)
      gphys_dse_ary.push( gphys_dse.rename("ml_dse_lat20_zonal").
			 set_att("line_name",$rezol).
			 set_att("ape_name", "dry static energy").
			 set_att("units","J kg-1") )
      gphys_mse_ary.push( gphys_mse.rename("ml_mse_lat20_zonal").
			 set_att("line_name",$rezol).
			 set_att("ape_name", "moist static energy").
			 set_att("units","J kg-1") )
      puts $rezol
    }
    mkfig_plot(gphys_dse_ary)
    mkfig_plot(gphys_mse_ary)

  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_tr_all_300day

    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_300day

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

      if item == "tr_tppn" 
	__gr_tr_tppn_300day
      elsif item == "tr_mslp"
	__gr_tr_mslp_300day
     elsif item == "U"
	__gr_tr_u250_300day
      elsif item == "V"
	__gr_tr_v250_300day
      elsif item == "OMG" 
	__gr_tr_om500_300day
      end
    }

  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.rename("tr_tppn_log").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].
    gphys = gphys[true,-401..-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").
    gphys = gphys[true,-401..-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").
    gphys = gphys[true,-401..-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").
    gphys = gphys[true,-401..-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").
#    gphys = gphys.stspct_fft("#{gphys.name}_log_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..80])
  end


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

  def __gr_tr_tppn_300day
    # 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,-801..-401].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 + 160).rename("time").
 put_att("units","days"))
    grid = Grid.new(grid_0,grid_1)
    gphys = GPhys.new(grid, gphys.data ).rename("tr_tppn_200")

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

    gphys = @data.gphys_open("tr_tppn").cut(true,0,true)
    gphys = gphys[true,-1201..-801].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 + 60).rename("time").
 put_att("units","days"))
    grid = Grid.new(grid_0,grid_1)
    gphys = GPhys.new(grid, gphys.data ).rename("tr_tppn_300")
    mkfig_plot(gphys.lon_lotate)
    # 時空間スペクトル
    __gr_spct(gphys)

  end

  def __gr_tr_mslp_300day
    # 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,-801..-401].set_lost_axes(lost_axis).
      rename("tr_mslp_200")
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate.
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)

    gphys = @data.gphys_open("tr_mslp").cut(true,0,true)
    gphys = gphys[true,-1201..-801].set_lost_axes(lost_axis). 
      rename("tr_mslp_300")
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate.
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)

  end

  def __gr_tr_u250_300day
    # 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,-801..-401].rename("tr_u250_200").
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate.
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)

    gphys = @data.gphys_open("U").cut(true,0.25,true)
    gphys = gphys[true,-1201..-801].rename("tr_u250_300").
      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_300day
    # 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,-801..-401].rename("tr_v250_200").
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate. 
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)

    gphys = @data.gphys_open("V").cut(true,0.25,true)
    gphys = gphys[true,-1201..-801].rename("tr_v250_300").
      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_300day
    # 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,-801..-401].rename("tr_om500_200").
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate. 
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)

    gphys = @data.gphys_open("OMG").cut(true,0.5,true)
    gphys = gphys[true,-1201..-801].rename("tr_om500_300").
      set_lost_axes(lost_axis)
    mkfig_plot((gphys - gphys.mean(0)).lon_lotate. 
	       add_lost_axes("(diff) from (mean) zonal"))
    # 時空間スペクトル
    __gr_spct(gphys)

  end

end