#!/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"  || item == "time_bnds" || item == "bounds_time" || item == "latitude" || item == "longitude" || item == "bounds_latitude" || item == "bounds_longitude" 

	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"|| $rezol == "CGAM" || $rezol == "DWD" || $rezol =~ /ECMWF/ || $rezol == "GFDL" || $rezol == "LASG" || $rezol == "MIT"|| $rezol == "MRI"
	  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" || item == "lat_bnds" || item == "lon_bnds" || item == "wt_lat_p" || item == "lat_p" || item == "lon_p" || item == "lat_p_bnds" || item == "lon_p_bnds" || item == "bnd" || item == "wt_lat" || item == "bounds_lat" || item == "bounds_lon"  || item == "longitude"  || item == "latitude" || item == "levelist" || item == "pres" || item == "time_bnds"  || item == "timebnd" || item == "lonbnd" || item == "latbnd" || item == "levbnd" || item == "pressure_bnds" || item == "latitude_bnds"  || item == "longitude_bnds" 


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

	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.rename("ml_phi") if  item == "ml_z" 

	if ( $rezol == "CSIRO" || $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "NCAR" || $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old" || $rezol == "MRI") && item == "ml_rh"
	  gphys = gphys*100
	end
	if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" || $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old"|| $rezol == "DWD" || $rezol == "MIT" || $rezol == "MRI" ) && 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[true,true,true,0] if $rezol =~ /ECMWF/

    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" ||  $rezol =~ /ECMWF/||  $rezol == "GFDL"
      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)
    if $rezol == "MRI"
      gphys_phi = @data.gphys_open("ml_z").mean(0)
    else
      gphys_phi = @data.gphys_open("ml_phi").mean(0)
    end

    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"  || $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old"|| $rezol == "DWD" || $rezol == "MIT" || $rezol == "MRI") 
      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}/org/#{$groupid}_TR_#{$expID}_a.nc",
			"#{$ncfile_path}/org/#{$groupid}_TR_#{$expID}_b.nc"])
      @data_tmp = Ape.new("#{$ncfile_path}/org/#{$groupid}_TR_#{$expID}_a.nc")
#      @data = Ape.new(["#{$ncfile_path}/org/ukmo_03f_TR_#{$expID}_a.nc",
#			"#{$ncfile_path}/org/ukmo_03f_TR_#{$expID}_b.nc"])
#      @data_tmp = Ape.new("#{$ncfile_path}/org/ukmo_03f_TR_#{$expID}_a.nc")
    elsif $rezol == "T159L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
    elsif $rezol == "T319L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc")
    elsif $rezol == "ECMWF_07"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR00_#{$expID}.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR00_#{$expID}.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|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
          
          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"))
          elsif $rezol == "T159L48_non" || $rezol == "T319L48_non"
            g = @data_tmp.gphys_open(name).cut(true,0,true)
            gphys = __T159L48_non_gphys_open(g,0)[true,-401..-1]
	  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"))
            grid_0 = gphys.grid_copy.axis(0)	
            grid = Grid.new(grid_0,grid_1)
            gphys = GPhys.new(grid, gphys.data ).set_lost_axes(gphys.lost_axes)
	  end

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

          if $rezol == "DWD" && gphys.name == "tr_mslp"
            gphys = gphys*100 
          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 - 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/

          gphys = @data.gphys_open(name)
  
	  # x-t ダイヤグラム
	  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"))
          elsif $rezol == "T159L48_non" || $rezol == "T319L48_non"
            g = @data_tmp.gphys_open(name).cut(true,0,true)
            gphys = __T159L48_non_gphys_open(g,0)[true,-401..-1]
	  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"))
            grid_0 = gphys.grid_copy.axis(0)	
            grid = Grid.new(grid_0,grid_1)
            gphys = GPhys.new(grid, gphys.data ).set_lost_axes(gphys.lost_axes)
	  end
	  
	  if $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")
    elsif $rezol == "T159L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
    elsif $rezol == "T319L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.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 ダイヤグラム
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).
            cut(true,-5..5,true)
          gphys = __T159L48_non_gphys_open(g,-5..5)[true,true,-401..-1].
            rename("#{name}_5deg_mean")
        else
          gphys = @data.gphys_open(name)
          gphys =  gphys.cut(true,-5..5,true)[true,true,-401..-1].
            rename("#{name}_5deg_mean")
        end

	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
        gphys = gphys*100 if $rezol =~ /ECMWF/ && gphys.name == "tr_mslp_5deg_mean"

	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 ダイヤグラム
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).
            cut(true,-5..5,true)
          gphys = __T159L48_non_gphys_open(g,-5..5)[true,true,-401..-1].
            rename("#{name}_5deg_mean")
        else
          gphys = @data.gphys_open(name)
          gphys =  gphys.cut(true,-5..5,true)[true,true,-401..-1].
            rename("#{name}_5deg_mean")
        end
	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 == "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")
    elsif $rezol == "T159L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
    elsif $rezol == "T319L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.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 ダイヤグラム
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).cut(true,0,true)
          gphys = __T159L48_non_gphys_open(g,0)[true,-801..-401].
            rename("#{name}_200")
        else
          gphys = @data.gphys_open(name)
          gphys =  gphys.cut(true,0,true)[true,-801..-401].
            rename("#{name}_200")
        end

	gphys = gphys.rename("tr_mslp_200")  if gphys.name == "tr_ps_200" 
        gphys = gphys*100 if $rezol =~ /ECMWF/ && gphys.name == "tr_mslp_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])

	#
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).cut(true,0,true)
          gphys = __T159L48_non_gphys_open(g,0)[true,-1201..-801].
            rename("#{name}_300")
        else
          gphys = @data.gphys_open(name)
          gphys =  gphys.cut(true,0,true)[true,-1201..-801].
            rename("#{name}_300")
        end

	gphys = gphys.rename("tr_mslp_300")  if gphys.name == "tr_ps_300" 
        gphys = gphys*100 if $rezol =~ /ECMWF/ && gphys.name == "tr_mslp_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 ダイヤグラム
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).cut(true,0,true)
          gphys = __T159L48_non_gphys_open(g,0)[true,-801..-401].
            rename("#{name}_200")
        else
          gphys = @data.gphys_open(name)
          gphys =  gphys.cut(true,0,true)[true,-801..-401].
            rename("#{name}_200")
        end
	if $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 ダイヤグラム
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).cut(true,0,true)
          gphys = __T159L48_non_gphys_open(g,0)[true,-1201..-801].
            rename("#{name}_300")
        else
          gphys = @data.gphys_open(name)
          gphys =  gphys.cut(true,0,true)[true,-1201..-801].
            rename("#{name}_300")
        end

	if $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")
    elsif $rezol == "T159L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc")
    elsif $rezol == "T319L48_non"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc")
      @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.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 ダイヤグラム
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).cut(true,0,true)
          gphys = __T159L48_non_gphys_open(g,0).
            rename("#{name}_360all")
        else
          gphys = @data.gphys_open(name)
          gphys =  gphys.cut(true,0,true).
            rename("#{name}_360all")
        end

	gphys = gphys.rename("tr_mslp_360all")  if gphys.name == "tr_ps_360all" 
        gphys = gphys*100 if $rezol =~ /ECMWF/ && gphys.name == "tr_mslp_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 ダイヤグラム
        if $rezol == "T159L48_non" || $rezol == "T319L48_non"
          g = @data_tmp.gphys_open(name).cut(true,0,true)
          gphys = __T159L48_non_gphys_open(g,0).
            rename("#{name}_360all")
        else
          gphys = @data.gphys_open(name)
          gphys = gphys.cut(true,0,true).
            rename("#{name}_360all")
        end
	if $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" || item == "bnd" || item == "wt_lat" || item == "wt_lon" || item == "lon_bnds" || item == "lat_bnds" || item == "wt_lat" || item == "wt_lat_u" || item == "lat_u"   || item == "lon_u"  || item == "lat_u_bnds"   || item == "lon_u_bnds" || item == "bounds_lat" || item == "bounds_lon"  || item == "longitude"  || item == "latitude"  || item == "gw" || item == "time_bnds" || item == "timebnd" || item == "lonbnd" || item == "latbnd"  || item == "latitude_73"  || item == "latitude_72" || item == "latitude_73_bnds" || item == "latitude_72_bnds"    || item == "bounds_time"  || item == "bounds_latitude" || item == "bounds_longitude"  || item == "bounds_latitude_73"  || item == "bounds_latitude_72" || item == "bounds_latitude_96" || item == "latitude_96" || item == "bounds_longitude_96" || item == "longitude_96" || item == "latitude_144" || item == "bounds_latitude_144" || item == "longitude_192" || item == "bounds_longitude_192"      

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

	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 == "CSIRO" || $rezol == "GSFC" || $rezol == "LASG" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" || $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old" || $rezol == "DWD" || $rezol =~ /ECMWF/ || $rezol == "GFDL" || $rezol == "MIT"|| $rezol == "MRI"
	  if item == "sh_tauv" || item == "sh_tauu"
	    gphys = -gphys
	  end
	end

	if $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "MGO" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" || $rezol == "CGAM"|| $rezol == "DWD" || $rezol =~ /ECMWF/ || $rezol == "GFDL"  || $rezol == "LASG" || $rezol == "MIT"|| $rezol == "MRI"
	  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[true,true,0] 	if $rezol =~ /ECMWF/

      gphys = gphys.lon_lotate
      mkfig_plot(gphys)

      # net ssw (ネットの短波放射)
      g1 = @data.gphys_open("sh_sswr")
      g2 = @data.gphys_open("sh_sswi")
      gphys = (g1 - g2).rename("sh_ssw").
	set_att("ape_name", "net_surface_shortwave_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[true,true,0] 	if $rezol =~ /ECMWF/

      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 == "nv" || item == "bnd" || item == "wt_lat" || item == "wt_lon" || item == "lon_bnds" || item == "lat_bnds" || item == "wt_lat" || item == "wt_lat_u" || item == "lat_u"   || item == "lon_u"  || item == "lat_u_bnds"   || item == "lon_u_bnds" || item == "bounds_lat" || item == "bounds_lon"  || item == "longitude"  || item == "latitude"  || item == "gw" || item == "time_bnds" || item == "timebnd" || item == "lonbnd" || item == "latbnd"  || item == "latitude_73"  || item == "latitude_72" || item == "latitude_73_bnds" || item == "latitude_72_bnds"    || item == "bounds_time"  || item == "bounds_latitude" || item == "bounds_longitude"  || item == "bounds_latitude_73"  || item == "bounds_latitude_72" || item == "bounds_latitude_96" || item == "latitude_96" || item == "bounds_longitude_96" || item == "longitude_96" || item == "latitude_144" || item == "bounds_latitude_144" || item == "longitude_192" || item == "bounds_longitude_192"      

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

	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" || $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old" || $rezol == "DWD" || $rezol =~ /ECMWF/ || $rezol == "GFDL" || $rezol == "MIT" || $rezol == "MRI"
	  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" || $rezol == "CGAM" || $rezol == "DWD"|| $rezol =~ /ECMWF/ || $rezol == "GFDL"|| $rezol == "LASG" || $rezol == "MIT" || $rezol == "MRI"
	  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"
      gphys = gphys[true,true,0] 	if $rezol =~ /ECMWF/
      mkfig_plot(gphys.mean(0))

      # net ssw (大気下端ネットの短波放射)
      g1 = @data.gphys_open("sh_sswi")
      g2 = @data.gphys_open("sh_sswr")
      gphys = (g1 - g2).rename("sh_ssw_zonal").
	set_att("ape_name", "net_surface_shortwave_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[true,true,0] 	if $rezol =~ /ECMWF/
      mkfig_plot(gphys.mean(0))


      # net ssw (ネットの短波放射=大気吸収量)
      g1 = @data.gphys_open("sh_sswr")
      g2 = @data.gphys_open("sh_sswi")
      g3 = @data.gphys_open("sh_sw_toai")
      g4 = @data.gphys_open("sh_sw_toar")
      gphys = (g1 - g2 + g3 - g4  ).rename("sh_sw_zonal").
	set_att("ape_name", "net_shortwave_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[true,true,0] 	if $rezol =~ /ECMWF/
      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

  
  #  #{$groupid}_PF_control.nc
  def nc_pf
    @data = Ape.new("#{$ncfile_path}#{$groupid}_PF_#{$expID}.nc")
    if $rezol =~ /ECMWF/ 
      @data = Ape.new("#{$ncfile_path}#{$groupid}_PFp_#{$expID}.nc")
    end

    @data.netcdf_open.var_names.each { |item|
      unless item == "lat" || item == "plev" || item == "lon" || item == "lev" || item == "sh_ps" || item == "ilev" || item == "pf_u_gwd" || item == "pf_v_gwd"|| item == "lat_bnds"|| item == "lon_bnds" || item == "hyai"|| item == "hybi" || item == "hyam" || item == "hybm" || item == "p0_ref" || item == "longitude" || item == "latitude" || item == "levelist" || item == "time"|| item == "pt"|| item == "gw" || item == "P0" || item == "time_bnds" || item == "timebnd" || item == "lonbnd"|| item == "latbnd" || item == "levbnd" || item == "lev_38" || item == "lev_38_bnds" || item == "latitude_73"  || item == "latitude_73_bnds" || item == "latitude_72" || item == "latitude_72_bnds" || item == "lev_12" || item == "lev_12_bnds" || item == "latitude_145"  || item == "latitude_145_bnds" || item == "latitude_144"  || item == "latitude_144_bnds" || item == "lev_38" || item == "lev_38_bnds"  || item == "longitude_bnds"

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

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

        if $rezol == "GSFC" && item == "pf_t_disp"
        elsif $rezol =~ /UKMO/ 
          if item == "pf_t_lw" || item == "pf_t_sw"
          else
            mkfig_plot(gphys)
          end
        else
          mkfig_plot(gphys.mean(0))
        end

      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" || item == "vertical" || item == "site"  || item == "time_bnds" || item == "bounds_time" || item == "latitude" || item == "longitude" || item == "bounds_latitude" || item == "bounds_longitude" 
# 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")


          p item
        if item == "gt_cppn" && $rezol =~ /non/
        elsif item == "gt_sw_toar" || item == "gt_sswr" || item == "gt_slwd" 
	  mkfig_plot(-gphys + gphys_con)
        else
          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" || item == "lat_bnds" || item == "lon_bnds" || item == "wt_lat_p" || item == "lat_p" || item == "lon_p" || item == "lat_p_bnds" || item == "lon_p_bnds" || item == "bnd" || item == "wt_lat" || item == "bounds_lat" || item == "bounds_lon"  || item == "longitude"  || item == "latitude" || item == "levelist" || item == "pres" || item == "time_bnds"  || item == "timebnd" || item == "lonbnd" || item == "latbnd" || item == "levbnd" || item == "pressure_bnds" || item == "latitude_bnds"  || item == "longitude_bnds" || item == "lev_bnds" # || item == "ml_rh" 

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

	gphys = gphys[true,true,true,0] if  $rezol =~ /ECMWF/
	gphys_con = gphys_con[true,true,true,0] if  $rezol =~ /ECMWF/

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

	gphys = gphys.rename("ml_phi") if  item == "ml_z" 
	gphys_con = gphys_con.rename("ml_phi") if item == "ml_z" 

	if ( $rezol == "CSIRO" || $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "NCAR" || $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old" || $rezol == "MRI" ) && item == "ml_rh"
	  gphys = gphys*100
	  gphys_con = gphys_con*100
	end
	if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" || $rezol == "CGAM"|| $rezol == "CSIRO_standard" || $rezol == "CSIRO_old" || $rezol == "DWD" || $rezol == "MIT" || $rezol == "MRI") && 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[true,true,true,0] if  $rezol =~ /ECMWF/

    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[true,true,true,0] if  $rezol =~ /ECMWF/

    gphys_con = gphys_con.mean(0).strm_function("ml_psi_anm").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1") 

    if  $rezol == "UKMO_n48" || $rezol == "UKMO_n96" ||  $rezol == "MGO" ||  $rezol =~ /ECMWF/ || $rezol == "GFDL"
      gphys = gphys*100
      gphys_con = gphys_con*100
    end

    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_t_con = @data_con.gphys_open("ml_t").mean(0)
    gphys_q_con = @data_con.gphys_open("ml_q").mean(0)
    if $rezol == "MRI"
      gphys_phi = @data.gphys_open("ml_z").mean(0)
      gphys_phi_con = @data_con.gphys_open("ml_z").mean(0)
    else
      gphys_phi = @data.gphys_open("ml_phi").mean(0)
      gphys_phi_con = @data_con.gphys_open("ml_phi").mean(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 =~ /ECMWF/ 
      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"|| $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old"|| $rezol == "DWD" || $rezol == "MIT"|| $rezol == "MRI") 
      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" || item == "bnd" || item == "wt_lat" || item == "wt_lon" || item == "lon_bnds" || item == "lat_bnds" || item == "wt_lat" || item == "wt_lat_u" || item == "lat_u"   || item == "lon_u"  || item == "lat_u_bnds"   || item == "lon_u_bnds" || item == "bounds_lat" || item == "bounds_lon"  || item == "longitude"  || item == "latitude"  || item == "gw" || item == "time_bnds" || item == "timebnd" || item == "lonbnd" || item == "latbnd"  || item == "latitude_73"  || item == "latitude_72" || item == "latitude_73_bnds" || item == "latitude_72_bnds"    || item == "bounds_time"  || item == "bounds_latitude" || item == "bounds_longitude"  || item == "bounds_latitude_73"  || item == "bounds_latitude_72" || item == "bounds_latitude_96" || item == "latitude_96" || item == "bounds_longitude_96" || item == "longitude_96" || item == "latitude_144" || item == "bounds_latitude_144" || item == "longitude_192" || item == "bounds_longitude_192"      

        p item
	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[true,true,0] if  $rezol =~ /ECMWF/
	gphys_con = gphys_con[true,true,0] if  $rezol =~ /ECMWF/


	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 == "CSIRO" || $rezol == "GSFC" || $rezol == "LASG" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" || $rezol == "CGAM" || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old" || $rezol == "DWD" || $rezol =~ /ECMWF/ || $rezol == "GFDL" || $rezol == "MIT"|| $rezol == "MRI"
	  if item == "sh_tauv" || item == "sh_tauu"
	    gphys = -gphys
	  end
	end

	if $rezol == "GSFC" || $rezol == "K1JAPAN" || $rezol == "MGO" || $rezol == "NCAR" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96" || $rezol == "CGAM" || $rezol == "DWD"  || $rezol =~ /ECMWF/   || $rezol == "GFDL" || $rezol == "LASG" || $rezol == "MIT"|| $rezol == "MRI"
	  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" || $rezol =~ /ECMWF/
	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_cut_t =  gphys.cut(true,0).rename("sh_tppn_lat0_anm"). 
	add_lost_axes(lost_axis)
      mkfig_plot(gphys_cut_t)

      gphys_mean_t =  gphys.cut(true,-30..30).mean(1).
        rename("sh_tppn_lat30m_anm"). 
	add_lost_axes(lost_axis)
      mkfig_plot(gphys_mean_t)

      gphys_mean_tt =  gphys.cut(true,-15..15).mean(1).
        rename("sh_tppn_lat15m_anm"). 
	add_lost_axes(lost_axis)
      mkfig_plot(gphys_mean_tt)

      # 赤道上 蒸発量 (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" || $rezol =~ /ECMWF/
	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_cut_e = gphys.cut(true,0).rename("sh_slh_lat0_anm").
	add_lost_axes(lost_axis)
      mkfig_plot(gphys_cut_e)

      gphys_mean_e = gphys.cut(true,-30..30).mean(1).rename("sh_slh_lat30m_anm").
	add_lost_axes(lost_axis)
      mkfig_plot(gphys_mean_e)

      gphys_mean_ee = gphys.cut(true,-15..15).mean(1).
        rename("sh_slh_lat15m_anm").
	add_lost_axes(lost_axis)
      mkfig_plot(gphys_mean_ee)

      gphys_mean_e = gphys_mean_e *0.4*10**(-6)
      gphys_mean_ee = gphys_mean_ee *0.4*10**(-6)
      gphys_cut_e = gphys_cut_e *0.4*10**(-6)
      
      gphys_mean_t = gphys_mean_t.rename("sh_tppn_evap_lat30m_anm").
	set_att("ape_name", "precipitation_and_evaporation_flux").
        set_att("line_name","sh_tppn")
      gphys_mean_e = gphys_mean_e.rename("sh_tppn_evap_lat30m_anm").
	set_att("ape_name", "precipitation_and_evaporation_flux").
        set_att("line_name","sh_evap")
      gphys_mean_tt = gphys_mean_tt.rename("sh_tppn_evap_lat15m_anm").
	set_att("ape_name", "precipitation_and_evaporation_flux").
        set_att("line_name","sh_tppn")
      gphys_mean_ee = gphys_mean_ee.rename("sh_tppn_evap_lat15m_anm").
	set_att("ape_name", "precipitation_and_evaporation_flux").
        set_att("line_name","sh_evap")
      gphys_cut_t = gphys_cut_t.rename("sh_tppn_evap_lat0_anm").
	set_att("ape_name", "precipitation_and_evaporation_flux").
        set_att("line_name","sh_tppn")
      gphys_cut_e = gphys_cut_e.rename("sh_tppn_evap_lat0_anm"). 
	set_att("ape_name", "precipitation_and_evaporation_flux").
        set_att("line_name","sh_evap")

      mkfig_plot([gphys_mean_t,gphys_mean_e])
      mkfig_plot([gphys_mean_tt,gphys_mean_ee])
      mkfig_plot([gphys_cut_t,gphys_cut_e])

    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 == "nv" || item == "bnd" || item == "wt_lat" || item == "wt_lon" || item == "lon_bnds" || item == "lat_bnds" || item == "wt_lat" || item == "wt_lat_u" || item == "lat_u"   || item == "lon_u"  || item == "lat_u_bnds"   || item == "lon_u_bnds" || item == "bounds_lat" || item == "bounds_lon"  || item == "longitude"  || item == "latitude"  || item == "gw" || item == "time_bnds" || item == "timebnd" || item == "lonbnd" || item == "latbnd"  || item == "latitude_73"  || item == "latitude_72" || item == "latitude_73_bnds" || item == "latitude_72_bnds"    || item == "bounds_time"  || item == "bounds_latitude" || item == "bounds_longitude"  || item == "bounds_latitude_73"  || item == "bounds_latitude_72" || item == "bounds_latitude_96" || item == "latitude_96" || item == "bounds_longitude_96" || item == "longitude_96" || item == "latitude_144" || item == "bounds_latitude_144" || item == "longitude_192" || item == "bounds_longitude_192"      
	
	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" || $rezol =~ /ECMWF/
	  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" || $rezol == "CGAM"|| $rezol == "CSIRO_standard" || $rezol == "CSIRO_old" || $rezol == "DWD" || $rezol =~ /ECMWF/   || $rezol == "GFDL" || $rezol == "MIT"|| $rezol == "MRI"
	  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" || $rezol == "CGAM" || $rezol == "DWD" || $rezol =~ /ECMWF/   || $rezol == "GFDL" || $rezol == "LASG" || $rezol == "MIT"|| $rezol == "MRI"
	  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"|| $rezol =~ /ECMWF/
	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" || $rezol =~ /ECMWF/ 
	gphys_ps = gphys_ps[true,true,0] 
	gphys_ps_con = gphys_ps_con[true,true,0] 
    end

    if  $rezol == "GFDL" || $rezol == "GSFC" 
      lost_axis = ["plev=92500 Pa",
        "(diff) from (mean) #{zonalcont}"]
    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,1000).lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,1000).lon_lotate

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

    elsif $rezol == "UKMO_n96" 
	gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	  cut(true,true,1000).lon_lotate
	gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	  cut(true,true,1000).lon_lotate

        gphys_u_con = @data_con.gphys_open("ml_u").
          cut(true,true,1000).lon_lotate
        gphys_v_con = @data_con.gphys_open("ml_v").
          cut(true,true,1000).lon_lotate
    elsif  $rezol =~ /ECMWF/ 
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,1000).lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,1000).lon_lotate

      gphys_u_con = @data_con.gphys_open("ml_u")[true,true,true,0].
        cut(true,true,1000).lon_lotate
      gphys_v_con = @data_con.gphys_open("ml_v")[true,true,true,0].
        cut(true,true,1000).lon_lotate
    elsif  $rezol == "GFDL" || $rezol == "GSFC" 
#      gphys_u = @data.gphys_open("ml_u").cut(true,true,1000).lon_lotate
#      gphys_v = @data.gphys_open("ml_v").cut(true,true,1000).lon_lotate
#      gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,1000).
#        lon_lotate
#      gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,1000).
#        lon_lotate
      gphys_u = @data.gphys_open("ml_u")[true,true,1].lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,1].lon_lotate
      gphys_u_con = @data_con.gphys_open("ml_u")[true,true,1].
        lon_lotate
      gphys_v_con = @data_con.gphys_open("ml_v")[true,true,1].
        lon_lotate

    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

      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
    end


    # これでいいのかなぁ, 欠損値処理

    if ( gphys_u.data.val.class == NArray )
      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)
    else 
      mask = gphys_u.gt(-900) * gphys_u.lt(10000)
      gphys_u = NArrayMiss.to_nam(gphys_u.data.val.to_na, mask.to_na)
      mask = gphys_u_con.gt(-900) * gphys_u_con.lt(10000)
      gphys_u_con = NArrayMiss.to_nam(gphys_u_con.data.val.to_na, mask.to_na)
      mask = gphys_v.gt(-900) * gphys_v.lt(10000)
      gphys_v = NArrayMiss.to_nam(gphys_v.data.val.to_na, mask.to_na)
      mask = gphys_v_con.gt(-900) * gphys_v_con.lt(10000)
      gphys_v_con = NArrayMiss.to_nam(gphys_v_con.data.val.to_na, mask.to_na)
    end


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

    if $rezol == "UKMO_n48" 
      if $expID == "1keq" 
        gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].
          cut(true,true,250)
        gphys_phi_con = @data_con.gphys_open("ml_phi").
          cut(true,true,250)
      else
        gphys_phi = @data.gphys_open("ml_phi").cut(true,true,250)
        gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,250)
      end
    elsif $rezol == "UKMO_n96" 	    
      if $expID == "control"
	gphys_phi = @data.gphys_open("ml_phi").cut(true,true,250)
      gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,250)
      else
	gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].
          cut(true,true,250)
      gphys_phi_con = @data_con.gphys_open("ml_phi"). 
          cut(true,true,250)
      end
    elsif $rezol =~ /ECMWF/ 	    
      gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].cut(true,true,250)
      gphys_phi_con = @data_con.gphys_open("ml_phi")[true,true,true,0].cut(true,true,250)
    elsif $rezol == "MRI"
      gphys_phi = @data.gphys_open("ml_z").cut(true,true,25000).rename("ml_phi")
      gphys_phi_con = @data_con.gphys_open("ml_z").cut(true,true,25000).rename("ml_phi")
    elsif $rezol == "GFDL"
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,250)
      gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,250)
    else
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,25000)
      gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,25000)
    end

    if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96"|| $rezol == "CGAM" || $rezol == "DWD"  || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old"  || $rezol == "MIT"|| $rezol == "MRI") 
      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" 
      if $expID == "1keq" 
        gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
          cut(true,true,250).lon_lotate
        gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
          cut(true,true,250).lon_lotate
        gphys_u_con = @data_con.gphys_open("ml_u").
          cut(true,true,250).lon_lotate
        gphys_v_con = @data_con.gphys_open("ml_v"). 
          cut(true,true,250).lon_lotate
      else
        gphys_u = @data.gphys_open("ml_u").
          cut(true,true,250).lon_lotate
        gphys_v = @data.gphys_open("ml_v").
          cut(true,true,250).lon_lotate
        gphys_u_con = @data_con.gphys_open("ml_u").
          cut(true,true,250).lon_lotate
        gphys_v_con = @data_con.gphys_open("ml_v").
          cut(true,true,250).lon_lotate
      end
    elsif $rezol == "UKMO_n96" 	    
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
        cut(true,true,250).lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
        cut(true,true,250).lon_lotate
      gphys_u_con = @data_con.gphys_open("ml_u").
        cut(true,true,250).lon_lotate
      gphys_v_con = @data_con.gphys_open("ml_v").
        cut(true,true,250).lon_lotate
    elsif $rezol =~ /ECMWF/ 	    
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,250).lon_lotate
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,250).lon_lotate
      gphys_u_con = @data_con.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,250).lon_lotate
      gphys_v_con = @data_con.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,250).lon_lotate
    elsif $rezol == "GFDL"
      gphys_u = @data.gphys_open("ml_u").cut(true,true,250).lon_lotate
      gphys_v = @data.gphys_open("ml_v").cut(true,true,250).lon_lotate

      gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,250).
        lon_lotate
      gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,250).
        lon_lotate
    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

      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
    end

    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
        if (gphys_u.data.val.class == NArray)
	  u_data[num,true]  =  gphys_u.data.val[num,true]
	  v_data[num,true]  =  gphys_v.data.val[num,true]
        else
	  u_data[num,true]  =  gphys_u.data.val.to_na[num,true]
	  v_data[num,true]  =  gphys_v.data.val.to_na[num,true]
        end
      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")

    if $rezol == "UKMO_n48" && $expID == "1keq" 
      gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].
	cut(true,true,850)
      gphys_phi_con = @data_con.gphys_open("ml_phi"). 
        cut(true,true,850)
    elsif $rezol == "UKMO_n96" 	    
      gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].
	cut(true,true,850)
      gphys_phi_con = @data_con.gphys_open("ml_phi"). 
        cut(true,true,850)
    elsif $rezol =~ /ECMWF/
      gphys_phi = @data.gphys_open("ml_phi")[true,true,true,0].
	cut(true,true,850)
      gphys_phi_con = @data_con.gphys_open("ml_phi")[true,true,true,0].
        cut(true,true,850)
    elsif $rezol == "GFDL"
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,850)
      gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,850)
    elsif $rezol == "MRI"
      gphys_phi = @data.gphys_open("ml_z").cut(true,true,85000).rename("ml_phi")
      gphys_phi_con = @data_con.gphys_open("ml_z").cut(true,true,85000).rename("ml_phi")
    else
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,85000)
      gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,85000)
    end

    if ( $rezol == "CSIRO" || $rezol == "MGO" || $rezol == " FRCGC" || $rezol == "UKMO_n48" || $rezol == "UKMO_n96"|| $rezol == "CGAM"|| $rezol == "DWD"  || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old"  || $rezol == "MIT"|| $rezol == "MRI") 
      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,850)
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,850)
      gphys_u_con = @data_con.gphys_open("ml_u").
	cut(true,true,850)
      gphys_v_con = @data_con.gphys_open("ml_v").
	cut(true,true,850)
    elsif $rezol == "UKMO_n96" 	    
	gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	  cut(true,true,850)
	gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	  cut(true,true,850)
        gphys_u_con = @data_con.gphys_open("ml_u").
          cut(true,true,850)
        gphys_v_con = @data_con.gphys_open("ml_v").
          cut(true,true,850)
    elsif $rezol =~ /ECMWF/ 	    
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
        cut(true,true,850)
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
        cut(true,true,850)
      gphys_u_con = @data_con.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,850)
      gphys_v_con = @data_con.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,850)
    elsif $rezol == "GFDL"
      gphys_u = @data.gphys_open("ml_u").cut(true,true,850)
      gphys_v = @data.gphys_open("ml_v").cut(true,true,850)
      gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,850)
      gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,850)
    else
      gphys_u = @data.gphys_open("ml_u").cut(true,true,85000)
      gphys_v = @data.gphys_open("ml_v").cut(true,true,85000)
      gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,85000)
      gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,85000)
    end

    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
        if (gphys_u.data.val.class == NArray)
	  u_data[num,true]  =  gphys_u.data.val[num,true]
	  v_data[num,true]  =  gphys_v.data.val[num,true]
        else 
	  u_data[num,true]  =  gphys_u.data.val.to_na[num,true]
	  v_data[num,true]  =  gphys_v.data.val.to_na[num,true]
        end 
      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,500)
      gphys_t_con = @data_con.gphys_open("ml_t").
	cut(true,true,500)
    elsif $rezol == "UKMO_n96" 	    
      gphys_t = @data.gphys_open("ml_t")[true,true,true,0].
	cut(true,true,500)
      gphys_t_con = @data_con.gphys_open("ml_t").
	cut(true,true,500)
    elsif $rezol =~ /ECMWF/
      gphys_t = @data.gphys_open("ml_t")[true,true,true,0].
	cut(true,true,500)
      gphys_t_con = @data_con.gphys_open("ml_t")[true,true,true,0].
	cut(true,true,500)
    elsif $rezol == "GFDL"
      gphys_t = @data.gphys_open("ml_t").cut(true,true,500)
      gphys_t_con = @data_con.gphys_open("ml_t").cut(true,true,500)
    else
      gphys_t = @data.gphys_open("ml_t").cut(true,true,50000)
      gphys_t_con = @data_con.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,500)
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,500)
      gphys_u_con = @data_con.gphys_open("ml_u").
	cut(true,true,500)
      gphys_v_con = @data_con.gphys_open("ml_v").
	cut(true,true,500)
    elsif $rezol == "UKMO_n96" 	    
	gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
	  cut(true,true,500)
	gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
	  cut(true,true,500)
        gphys_u_con = @data_con.gphys_open("ml_u").
          cut(true,true,500)
        gphys_v_con = @data_con.gphys_open("ml_v").
          cut(true,true,500)
    elsif $rezol =~ /ECMWF/
      gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
        cut(true,true,500)
      gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
        cut(true,true,500)
      gphys_u_con = @data_con.gphys_open("ml_u")[true,true,true,0].
	cut(true,true,500)
      gphys_v_con = @data_con.gphys_open("ml_v")[true,true,true,0].
	cut(true,true,500)
    elsif $rezol == "GFDL"
      gphys_u = @data.gphys_open("ml_u").cut(true,true,500)
      gphys_v = @data.gphys_open("ml_v").cut(true,true,500)
      gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,500)
      gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,500)
    else
      gphys_u = @data.gphys_open("ml_u").cut(true,true,50000)
      gphys_v = @data.gphys_open("ml_v").cut(true,true,50000)

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

    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
        if (gphys_u.data.val.class == NArray)
	  u_data[num,true]  =  gphys_u.data.val[num,true]
	  v_data[num,true]  =  gphys_v.data.val[num,true]
        else 
	  u_data[num,true]  =  gphys_u.data.val.to_na[num,true]
	  v_data[num,true]  =  gphys_v.data.val.to_na[num,true]
        end
      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")
    gphys_t_con = @data_con.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

    if $rezol =~ /ECMWF/ 	    
      gphys_t = gphys_t[true,true,true,0] 
      gphys_t_con = gphys_t_con[true,true,true,0] 
    end

    gphys_t = gphys_t.cut(true,0,true)
    gphys_t_con = gphys_t_con.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")
    gphys_u_con = @data_con.gphys_open("ml_u")
    gphys_om_con = @data_con.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

    if $rezol =~ /ECMWF/ 	    
      gphys_u = gphys_u[true,true,true,0] 
      gphys_om = gphys_om[true,true,true,0] 
      gphys_u_con = gphys_u_con[true,true,true,0] 
      gphys_om_con = gphys_om_con[true,true,true,0] 
    end

    gphys_u = gphys_u.cut(true,0,true)
    gphys_om = gphys_om.cut(true,0,true)

    gphys_u_con = gphys_u_con.cut(true,0,true)
    gphys_om_con = gphys_om_con.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
        if ( gphys_u.data.val.class == NArray ) 
  	  u_data[num,true]  =  gphys_u.data.val[num,true]
	  om_data[num,true]  = gphys_om.data.val[num,true]
        else
	  u_data[num,true]  =  gphys_u.data.val.to_na[num,true]
	  om_data[num,true]  = gphys_om.data.val.to_na[num,true]
        end
      elsif   
	u_data[num,true]  = 0.0
	om_data[num,true]  = 0.0
      end
    }  
    
    if $rezol == "GFDL" 	    
      u_data[true,0]  = 0.0
      om_data[true,0] = 0.0
    end

    mkfig_plot(gphys_t,u_data,-om_data)
    
  end

  # Takaya Nakamura flux ---------------------------------------------
  
  def nc_sh_tnflux_level(zonalcont="control")

    lost_axis = ["TNflux: (basic state) #{zonalcont}"]

    dir = "/home/yukiko/work/ape/NetCDF/other_group/tnfx_netcdf/"

    [925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30].
      each{ |level|

      @data = Ape.new("#{dir}/TNX_#{$rezol}_#{$expID}.ctl")
      gphys_u = @data.gphys_open("TNX")[true,true,true,0].
        cut(true,true,level).lon_lotate
      @data = Ape.new("#{dir}/TNY_#{$rezol}_#{$expID}.ctl")
      gphys_v = @data.gphys_open("TNY")[true,true,true,0].
        cut(true,true,level).lon_lotate
      @data = Ape.new("#{dir}/TNZ_#{$rezol}_#{$expID}.ctl")
      gphys_z = @data.gphys_open("TNZ")[true,true,true,0].
        cut(true,true,level).lon_lotate
      @data = Ape.new("#{dir}/DIV_#{$rezol}_#{$expID}.ctl")
      gphys_div = @data.gphys_open("DIV")[true,true,true,0].
        cut(true,true,level).lon_lotate

      gphys_z = gphys_z.
        set_lost_axes(lost_axis).
        add_lost_axes("z=#{level} hPa").
        rename("sh_tnztn#{level}").
        set_att("ape_name","TN flux z, (x,y)").
        set_att("units","m2 s-2, (m2 s-2, m2 s-2)")
      
      gphys_div = gphys_div.
        set_lost_axes(lost_axis).
        add_lost_axes("z=#{level} hPa").
        rename("sh_tndivtn#{level}").
        set_att("ape_name","TN flux div, (x,y)").
        set_att("units","m day-1, (m2 s-2, m2 s-2)")
      
      @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")
      @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      if $rezol == "GFDL"
        gphys_phi = @data.gphys_open("ml_phi").cut(true,true,level)
        gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,level)
      else
        gphys_phi = @data.gphys_open("ml_phi").cut(true,true,level*100)
        gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,level*100)
      end
      
      lost_axis_2 = ["phi: (diff) from (mean) #{zonalcont}"]
      gphys_phi = ( gphys_phi - gphys_phi_con.mean(0) ).lon_lotate.
        set_lost_axes(lost_axis).add_lost_axes("z=#{level} hPa").
        add_lost_axes(lost_axis_2).
        rename("sh_phitn#{level}").
        set_att("ape_name","phi, TN flux (x,y)").
        set_att("units","m, (m2 s-2, m2 s-2)")

      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)
      
      u_data_tmp  = NArray.sfloat( gphys_u.axis(0).pos.val.size ,
                                   gphys_u.axis(1).pos.val.size ).fill!(0.0)
      v_data_tmp  = 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
          if (gphys_u.data.val.class == NArray)
            u_data[num,true]  =  gphys_u.data.val[num,true]
            v_data[num,true]  =  gphys_v.data.val[num,true]
            u_data_tmp[num,true]  =  gphys_u.data.val[num,true]
            v_data_tmp[num,true]  =  gphys_v.data.val[num,true]
          else
            u_data[num,true]  =  gphys_u.data.val.to_na[num,true]
            v_data[num,true]  =  gphys_v.data.val.to_na[num,true]
            u_data_tmp[num,true]  =  gphys_u.data.val.to_na[num,true]
            v_data_tmp[num,true]  =  gphys_v.data.val.to_na[num,true]
          end
        elsif   
          u_data[num,true]  = 0.0
          v_data[num,true]  = 0.0
          u_data_tmp[num,true]  = 0.0
          v_data_tmp[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
      }
      
      miss = 10
      gphys_u.axis(0).pos.val.size.times{ |numu|
        gphys_u.axis(1).pos.val.size.times{ |numv|
          if (u_data_tmp[numu,numv] < -miss) || (u_data_tmp[numu,numv] > miss) || (v_data_tmp[numu,numv] < -miss) || (v_data_tmp[numu,numv] > miss )
            u_data[numu,numv] = 0.0
            v_data[numu,numv] = 0.0
#          u_data_tmp[numu,numv] = 0.0
#          v_data_tmp[numu,numv] = 0.0
          end
        }
      }
      mkfig_plot(gphys_phi,u_data,v_data)
      mkfig_plot(gphys_z,u_data,v_data)
      mkfig_plot(gphys_div,u_data,v_data)
      
      size= (gphys_u.axis(1).pos.val.size/ 2).to_i
      gphys_phi = gphys_phi[true,size..-1].rename("sh_phitn#{level}_lath")
    gphys_div = gphys_div[true,size..-1].rename("sh_tndivtn#{level}_lath")
    gphys_z  = gphys_z[true,size..-1].rename("sh_tnztn#{level}_lath")

    u_data = u_data_tmp[true,size..-1]
    v_data = v_data_tmp[true,size..-1]

    u_data[true,0..2] = 0.0
    v_data[true,0..2] = 0.0

    u_data[true,-2..-1] = 0.0
    v_data[true,-2..-1] = 0.0

    miss = 20
    gphys_u.axis(0).pos.val.size.times{ |numu|
      gphys_u.axis(1).pos.val.size.times{ |numv|
        if (u_data_tmp[numu,numv] < -miss) || (u_data_tmp[numu,numv] > miss) || (v_data_tmp[numu,numv] < -miss) || (v_data_tmp[numu,numv] > miss )
          u_data_tmp[numu,numv] = 0.0
          v_data_tmp[numu,numv] = 0.0
        end
      }
    }

    mkfig_plot(gphys_phi,u_data,v_data)
    mkfig_plot(gphys_z,u_data,v_data)
    mkfig_plot(gphys_div,u_data,v_data)

}
    
  end
  
  def nc_sh_tnflux(zonalcont="control")

    lost_axis = ["TNflux: (basic state) #{zonalcont}"]

    dir = "/home/yukiko/work/ape/NetCDF/other_group/tnfx_netcdf/"
    @data = Ape.new("#{dir}/TNX_#{$rezol}_#{$expID}.ctl")
    gphys_u = @data.gphys_open("TNX")[true,true,true,0].
      cut(true,true,250).lon_lotate
    @data = Ape.new("#{dir}/TNY_#{$rezol}_#{$expID}.ctl")
    gphys_v = @data.gphys_open("TNY")[true,true,true,0].
      cut(true,true,250).lon_lotate
    @data = Ape.new("#{dir}/TNZ_#{$rezol}_#{$expID}.ctl")
    gphys_z = @data.gphys_open("TNZ")[true,true,true,0].
      cut(true,true,850).lon_lotate
    @data = Ape.new("#{dir}/DIV_#{$rezol}_#{$expID}.ctl")
    gphys_div = @data.gphys_open("DIV")[true,true,true,0].
      cut(true,true,250).lon_lotate

    gphys_z = gphys_z.
      set_lost_axes(lost_axis).
      add_lost_axes("x,y: z=250 hPa, z: z=850 hPa").
      rename("sh_tnz850tn250").
      set_att("ape_name","TN flux z, (x,y)").
      set_att("units","m2 s-2, (m2 s-2, m2 s-2)")

    gphys_div = gphys_div.
      set_lost_axes(lost_axis).
      add_lost_axes("z=250 hPa").
      rename("sh_tndiv250tn250").
      set_att("ape_name","TN flux div, (x,y)").
      set_att("units","m day-1, (m2 s-2, m2 s-2)")

    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{zonalcont}.nc")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
    if $rezol == "GFDL"
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,250)
      gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,250)
    else
      gphys_phi = @data.gphys_open("ml_phi").cut(true,true,25000)
      gphys_phi_con = @data_con.gphys_open("ml_phi").cut(true,true,25000)
    end

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

    lost_axis_2 = ["phi: (diff) from (mean) #{zonalcont}"]
    gphys_phi = ( gphys_phi - gphys_phi_con.mean(0) ).lon_lotate.
      set_lost_axes(lost_axis).add_lost_axes("z=250 hPa").
      add_lost_axes(lost_axis_2).
      rename("sh_phitn250").
      set_att("ape_name","phi, TN flux (x,y)").
      set_att("units","m, (m2 s-2, m2 s-2)")

    lost_axis_2 = ["ps: (diff) from (mean) #{zonalcont}"]
    gphys_ps = ( gphys_ps - gphys_ps_con.mean(0) ).lon_lotate.
      set_lost_axes(lost_axis).
      add_lost_axes(lost_axis_2).
      rename("sh_pstn250").
      set_att("ape_name","ps, TN flux(x,y)").
      set_att("units","Pa, (m2 s-2, m2 s-2)")

    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)

    u_data_tmp  = NArray.sfloat( gphys_u.axis(0).pos.val.size ,
			    gphys_u.axis(1).pos.val.size ).fill!(0.0)
    v_data_tmp  = 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
        if (gphys_u.data.val.class == NArray)
	  u_data[num,true]  =  gphys_u.data.val[num,true]
	  v_data[num,true]  =  gphys_v.data.val[num,true]
	  u_data_tmp[num,true]  =  gphys_u.data.val[num,true]
	  v_data_tmp[num,true]  =  gphys_v.data.val[num,true]
        else
	  u_data[num,true]  =  gphys_u.data.val.to_na[num,true]
	  v_data[num,true]  =  gphys_v.data.val.to_na[num,true]
	  u_data_tmp[num,true]  =  gphys_u.data.val.to_na[num,true]
	  v_data_tmp[num,true]  =  gphys_v.data.val.to_na[num,true]
        end
      elsif   
	u_data[num,true]  = 0.0
	v_data[num,true]  = 0.0
	u_data_tmp[num,true]  = 0.0
	v_data_tmp[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
    }

    miss = 10
    gphys_u.axis(0).pos.val.size.times{ |numu|
      gphys_u.axis(1).pos.val.size.times{ |numv|
        if (u_data_tmp[numu,numv] < -miss) || (u_data_tmp[numu,numv] > miss) || (v_data_tmp[numu,numv] < -miss) || (v_data_tmp[numu,numv] > miss )
          u_data[numu,numv] = 0.0
          v_data[numu,numv] = 0.0
#          u_data_tmp[numu,numv] = 0.0
#          v_data_tmp[numu,numv] = 0.0
        end
      }
    }
    mkfig_plot(gphys_phi,u_data,v_data)
    mkfig_plot(gphys_ps,u_data,v_data)
    mkfig_plot(gphys_z,u_data,v_data)
    mkfig_plot(gphys_div,u_data,v_data)

    size= (gphys_u.axis(1).pos.val.size/ 2).to_i
    gphys_phi = gphys_phi[true,size..-1].rename("sh_phitn250_lath")
    gphys_ps = gphys_ps[true,size..-1].rename("sh_pstn250_lath")
    gphys_div = gphys_div[true,size..-1].rename("sh_tndiv250tn250_lath")
    gphys_z  = gphys_z[true,size..-1].rename("sh_tnz850tn250_lath")

    u_data = u_data_tmp[true,size..-1]
    v_data = v_data_tmp[true,size..-1]

    u_data[true,0..2] = 0.0
    v_data[true,0..2] = 0.0

    u_data[true,-2..-1] = 0.0
    v_data[true,-2..-1] = 0.0

    miss = 20
    gphys_u.axis(0).pos.val.size.times{ |numu|
      gphys_u.axis(1).pos.val.size.times{ |numv|
        if (u_data_tmp[numu,numv] < -miss) || (u_data_tmp[numu,numv] > miss) || (v_data_tmp[numu,numv] < -miss) || (v_data_tmp[numu,numv] > miss )
          u_data_tmp[numu,numv] = 0.0
          v_data_tmp[numu,numv] = 0.0
        end
      }
    }

    mkfig_plot(gphys_phi,u_data,v_data)
    mkfig_plot(gphys_ps,u_data,v_data)
    mkfig_plot(gphys_z,u_data,v_data)
    mkfig_plot(gphys_div,u_data,v_data)

    
  end
  
  # 渦度 -----------------------------------------------------

  # #{$groupid}_ML_control.nc
  def nc_sh_vor_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")

    [1000,850,500,250].each{ |plev|
      
      if $rezol == "UKMO_n48" 
        if $expID == "1keq" 
          gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
            cut(true,true,plev).lon_lotate
          gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
            cut(true,true,plev).lon_lotate
          gphys_u_con = @data_con.gphys_open("ml_u").
            cut(true,true,plev).lon_lotate
          gphys_v_con = @data_con.gphys_open("ml_v"). 
            cut(true,true,plev).lon_lotate
        else
          gphys_u = @data.gphys_open("ml_u").
            cut(true,true,plev).lon_lotate
          gphys_v = @data.gphys_open("ml_v").
            cut(true,true,plev).lon_lotate
          gphys_u_con = @data_con.gphys_open("ml_u").
            cut(true,true,plev).lon_lotate
          gphys_v_con = @data_con.gphys_open("ml_v").
            cut(true,true,plev).lon_lotate
        end
      elsif $rezol == "UKMO_n96" 	    
        gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
          cut(true,true,plev).lon_lotate
        gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
          cut(true,true,plev).lon_lotate
        gphys_u_con = @data_con.gphys_open("ml_u").
          cut(true,true,plev).lon_lotate
        gphys_v_con = @data_con.gphys_open("ml_v").
          cut(true,true,plev).lon_lotate
      elsif $rezol =~ /ECMWF/ 	    
        gphys_u = @data.gphys_open("ml_u")[true,true,true,0].
          cut(true,true,plev).lon_lotate
        gphys_v = @data.gphys_open("ml_v")[true,true,true,0].
          cut(true,true,plev).lon_lotate
        gphys_u_con = @data_con.gphys_open("ml_u")[true,true,true,0].
          cut(true,true,plev).lon_lotate
        gphys_v_con = @data_con.gphys_open("ml_v")[true,true,true,0].
          cut(true,true,plev).lon_lotate
      elsif $rezol == "GFDL"
        gphys_u = @data.gphys_open("ml_u").cut(true,true,plev).lon_lotate
        gphys_v = @data.gphys_open("ml_v").cut(true,true,plev).lon_lotate
        
        gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,plev).
          lon_lotate
        gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,plev).
         lon_lotate
      else
        gphys_u = @data.gphys_open("ml_u").cut(true,true,plev*100).lon_lotate
        gphys_v = @data.gphys_open("ml_v").cut(true,true,plev*100).lon_lotate
        
        gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,plev*100).
          lon_lotate
        gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,plev*100).
        lon_lotate
      end

      gphys = gphys_u.copy
      gphys_con = gphys_u.copy

      lat = gphys_u.grid_copy.axis(1)

      lat.pos.val.size.times{ |num| 

        lon = gphys_u.grid_copy.axis(0)
        lat = gphys_u.grid_copy.axis(1)
        
        lon =  2.0 * 3.14159 * 6371.0 * 1000.0 * lon.pos.val / 360  * cos(2.0*3.14159*lat.pos.val[num]/360.0)
        lat =  2.0 * 3.14159 * 6371.0 * 1000.0 * lat.pos.val / 360

        lon = Axis.new().set_pos(VArray.new(lon).rename("lon").put_att("units","degrees_east")) 
        lat = Axis.new().set_pos(VArray.new(lat).rename("lat").put_att("units","degrees_north")) 

        grid = Grid.new(lon, lat)
        gu = GPhys.new( grid, gphys_u.data )
        gv = GPhys.new( grid, gphys_v.data )

        guc = GPhys.new( grid, gphys_u_con.data )
        gvc = GPhys.new( grid, gphys_v_con.data )

        gphys[true,num] = 
          ( GPhys::Derivative::threepoint_O2nd_deriv(gv, 0) -
              GPhys::Derivative::threepoint_O2nd_deriv(gu, 1) )[true,num].
          data.val
        
        gphys_con[true,num] =
          ( GPhys::Derivative::threepoint_O2nd_deriv(gvc, 0) -
              GPhys::Derivative::threepoint_O2nd_deriv(guc, 1) )[true,num].
          data.val

      }

      gphys = ( gphys-gphys_con.mean(0) ).
        rename("sh_vor#{plev}_anm").
        set_att("ape_name","vorticity").
        set_att("units","s-1").
        add_lost_axes(lost_axis)

      mkfig_plot( gphys )

    }

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

            if idsize == 0
              gphys = @data.gphys_open(item).
                rename("#{item}").
                set_att("line_name",$expID)
            else 
              gphys = @data.gphys_open(item).
                rename("#{item}_#{idsize}").
                set_att("line_name",$expID)
            end
	    gphys_ary.push( gphys )
	  }
	  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")
        if idsize == 0
          gphys = gphys.
            rename("#{gphys.name}").
            set_att("line_name",$expID) 
        else 
          gphys = gphys.
            rename("#{gphys.name}_#{idsize}").
            set_att("line_name",$expID)
        end
        gphys_ary.push( gphys )
      }
      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( id )
    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
            if idsize == 0
              gphys_ary.push(gphys.
                               rename("#{gphys.name}_zonal").
                               set_att("line_name",$expID) )
            else
              gphys_ary.push(gphys.
                               rename("#{gphys.name}_zonal_#{idsize}").
                               set_att("line_name",$expID) )
            end
	  }
	  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
        if idsize == 0
          gphys_ary.push(gphys.mean(0).
                           rename("#{gphys.name}").
                           set_att("line_name",$expID) )
        else
          gphys_ary.push(gphys.mean(0).
                           rename("#{gphys.name}_#{idsize}").
                           set_att("line_name",$expID) )
        end
      }

      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}_ML_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 = "/home/yukiko/work/ape/NetCDF/#{$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",$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 = "/home/yukiko/work/ape/NetCDF/#{$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"|| $rezol == "CGAM"|| $rezol == "DWD"  || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old"  || $rezol == "MIT"|| $rezol == "MRI"  ) 
	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}"

    if $expID == "control" 
      rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}","T319L48_#{$cumulus}"]
    else
      rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}"]
    end

    $ncfile_path = "/home/yukiko/work/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 = "/home/yukiko/work/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 = "/home/yukiko/work/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_rezol
    $file_label = "#{$cumulus}_#{$expID}"

    if $expID == "control" && $cumulus == "eml"
      rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}","T319L48_#{$cumulus}"]
    else
    rezol = ["T39L48_#{$cumulus}", "T39L24_#{$cumulus}", "T39L96_#{$cumulus}", "T79L48_#{$cumulus}", "T159L48_#{$cumulus}"]
    end

    $ncfile_path = "/home/yukiko/work/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 = "/home/yukiko/work/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 = "/home/yukiko/work/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_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/work/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 = "/home/yukiko/work/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" || $rezol == "CGAM"  || $rezol == "CSIRO_standard" || $rezol == "CSIRO_old"|| $rezol == "DWD" || $rezol == "MIT" || $rezol == "MRI") 
	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 = "/home/yukiko/work/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 = "/home/yukiko/work/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 = "/home/yukiko/work/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 = "/home/yukiko/work/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 = "/home/yukiko/work/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 = "/home/yukiko/work/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}_ML_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 = "/home/yukiko/work/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 = "/home/yukiko/work/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"|| $rezol == "CGAM" || $rezol == "DWD" || $rezol == "CSIRO_old"  || $rezol == "CSIRO_standard"  || $rezol == "MIT"|| $rezol == "MRI" ) 
	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 __T159L48_non_gphys_open(g,dimcut)

    grid0size = g.grid_copy.axis(0).pos.val.size
    grid_0    = g.grid_copy.axis(0)

    unless dimcut == 0
      grid1size = g.grid_copy.axis(1).pos.val.size
      grid_1    = g.grid_copy.axis(1)
    end

    timesize = 4*360
    grid_time =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(timesize).indgen!/4+ 0.25).rename("time").
    put_att("units","days"))

    if $rezol == "T159L48_non"
      if dimcut == 0
        data =  NArray.sfloat(grid0size,timesize).fill!(0.0)
        grid = Grid.new(grid_0, grid_time)

        data[true,0..(2*360-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_half1.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(2*360)..-1] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc",
                            g.name).cut(true,dimcut,true).data.val
      else
        data =  NArray.sfloat(grid0size,grid1size,timesize).fill!(0.0)
        grid = Grid.new(grid_0, grid_1, grid_time)

        data[true,true,0..(2*360-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_half1.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(2*360)..-1] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_half2.nc",
                            g.name).cut(true,dimcut,true).data.val
      end

    elsif $rezol == "T319L48_non"

      if dimcut == 0
        data =  NArray.sfloat(grid0size,timesize).fill!(0.0)
        grid = Grid.new(grid_0, grid_time)

        data[true,0..(4*30-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0031.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30)..(4*30*2-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0032.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*2)..(4*30*3-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0033.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*3)..(4*30*4-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0034.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*4)..(4*30*5-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0035.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*5)..(4*30*6-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0036.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*6)..(4*30*7-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0037.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*7)..(4*30*8-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0038.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*8)..(4*30*9-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0039.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*9)..(4*30*10-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0040.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*10)..(4*30*11-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0041.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,(4*30*11)..(4*30*12-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc",
                            g.name).cut(true,dimcut,true).data.val
      else
        data =  NArray.sfloat(grid0size,grid1size,timesize).fill!(0.0)
        grid = Grid.new(grid_0, grid_1, grid_time)

        data[true,true,0..(4*30-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0031.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30)..(4*30*2-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0032.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*2)..(4*30*3-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0033.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*3)..(4*30*4-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0034.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*4)..(4*30*5-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0035.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*5)..(4*30*6-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0036.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*6)..(4*30*7-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0037.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*7)..(4*30*8-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0038.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*8)..(4*30*9-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0039.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*9)..(4*30*10-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0040.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*10)..(4*30*11-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0041.nc",
                            g.name).cut(true,dimcut,true).data.val
        data[true,true,(4*30*11)..(4*30*12-1)] =
          GPhys::NetCDF_IO.open(
                        "#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc",
                            g.name).cut(true,dimcut,true).data.val
      end
    end

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

    return gphys

  end

  def __gr2_gphys_make(g,dimcut)

    grid0size = g.grid_copy.axis(0).pos.val.size
    grid_0    = g.grid_copy.axis(0)
    timesize = 360*4
    grid_time =
        Axis.new().
        set_pos(VArray.new(NArray.sfloat(timesize).indgen!/4+ 0.25).rename("time").
      put_att("units","days"))

    data =  NArray.sfloat(grid0size,timesize).fill!(0.0)
    grid = Grid.new(grid_0, grid_time)

    if $rezol =~ /T159L48/

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

    elsif $rezol =~ /T319L48_eml/

      data[true,(-30*4*4)..(-30*4*3-1)] =
        GPhys::NetCDF_IO.open(
                            "#{$gr2ncfile_path}#{$rezol}_expID01_0039.nc",
                            g.name).cut(true,dimcut,true).data.val
      data[true,(-30*4*3)..(-30*4*2-1)] =
        GPhys::NetCDF_IO.open(
                            "#{$gr2ncfile_path}#{$rezol}_expID01_0040.nc",
                            g.name).cut(true,dimcut,true).data.val
      data[true,(-30*4*2)..(-30*4-1)] =
        GPhys::NetCDF_IO.open(
                            "#{$gr2ncfile_path}#{$rezol}_expID01_0041.nc",
                            g.name).cut(true,dimcut,true).data.val
      data[true,(-30*4)..-1] =
        GPhys::NetCDF_IO.open(
                            "#{$gr2ncfile_path}#{$rezol}_expID01_0042.nc",
                            g.name).cut(true,dimcut,true).data.val

    else

      data[true,0..(6*30*4-1)] =
        GPhys::NetCDF_IO.open(
                            "#{$gr2ncf