#!/usr/bin/env ruby



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

class Ape_mkfig_agcm

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

  def Ape_mkfig_agcm.help
    print  <<EOF

  ==========================================================
   class Ape_mkfig_agcm のメソッド (lib-ape-agcm5-view.rb) 
  ---------------------------------------------------------
   * Ape_mkfig_agcm.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_agcm.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
=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

  end

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

  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
    }
  end

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

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

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

    # 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
    @data = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |name|
      if name =~ /lw_toa|mslp|500|250/
	if name =~ /lw_toa|mslp/
	  lost_axis = []
	elsif name =~ /250/
	  lost_axis = ["plev=25000 Pa"]
	elsif name =~ /500/
	  lost_axis = ["plev=50000 Pa"]
	end

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

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

      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
	gphys = gphys.cut(true,0,true)[true,-400..-1]
	gphys = gphys.lon_lotate unless $expID =~ /H1998/
	mkfig_plot(gphys)

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

  end

  #  #{$groupid}_SH_control.nc
  def nc_sh
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item == "lon" || item == "lat" 
	gphys = @data.gphys_open(item)
	gphys = gphys.lon_lotate unless $expID =~ /H1998/
	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" 
	gphys = @data.gphys_open(item)
	gphys = gphys.rename("#{gphys.name}_zonal")
	mkfig_plot(gphys.mean(0))
      end
    }

  end
  
  #  #{$groupid}_MF_control.nc
  def nc_mf
    @data = Ape.new("#{$ncfile_path}#{$groupid}_MF_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item == "lat" || item == "plev" 
	gphys = @data.gphys_open(item)
	mkfig_plot(gphys)
      end
    }
  end

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

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

  end

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

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

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

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

    # 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" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) zonal of #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	name = gphys.name
	gphys = (gphys - gphys_con.mean(0)).rename("#{gphys.name}_anm")
	gphys = gphys.lon_lotate unless $expID =~ /H1998/
	mkfig_plot(gphys)

	# 赤道上の値
	gphys = gphys.set_lost_axes("")
	gphys =  gphys.cut(true,0).rename("#{name}_lat0_anm"). 
	  add_lost_axes(lost_axis)
	mkfig_plot(gphys)
      end
    }

    # 赤道上の値 (tppn,evap)
    lost_axis = ["(diff) from (mean) zonal of #{zonalcont}"]
    gphys_con = @data_con.gphys_open("sh_tppn")
    gphys = @data.gphys_open("sh_tppn")
    gphys_t = (gphys - gphys_con.mean(0))
    gphys_t = gphys_t.set_lost_axes("")
    gphys_cut_t =  gphys_t.cut(true,0).rename("sh_tppn_evap_lat0_anm"). 
      add_lost_axes(lost_axis)
    gphys_mean_t =  gphys_t.cut(true,-30..30).mean(1).
      rename("sh_tppn_evap_lat30m_anm"). 
      add_lost_axes(lost_axis)
    gphys_mean_tt =  gphys_t.cut(true,-15..15).mean(1).
      rename("sh_tppn_evap_lat15m_anm"). 
      add_lost_axes(lost_axis)
    
    gphys_con = @data_con.gphys_open("sh_evap")
    gphys = @data.gphys_open("sh_evap")
    gphys_e = (gphys - gphys_con.mean(0))
    gphys_e = gphys_e.set_lost_axes("")
    gphys_cut_e =  gphys_e.cut(true,0).rename("sh_tppn_evap_lat0_anm")
    gphys_mean_e =  gphys_e.cut(true,-30..30).mean(1).
      rename("sh_tppn_evap_lat30m_anm"). 
      add_lost_axes(lost_axis)
    gphys_mean_ee =  gphys_e.cut(true,-15..15).mean(1).
      rename("sh_tppn_evap_lat15m_anm"). 
      add_lost_axes(lost_axis)

    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_cut_t,gphys_cut_e])
    mkfig_plot([gphys_mean_t,gphys_mean_e])
    mkfig_plot([gphys_mean_tt,gphys_mean_ee])

  end

  def nc_sh_zonal_anm(zonalcont="control")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{zonalcont}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item == "lon" || item == "lat" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_zonal_anm")
	mkfig_plot(gphys.mean(0)- gphys_con.mean(0))
      end
    }

  end
  
  def nc_mf_anm(zonalcont="control")
    @data = Ape.new("#{$ncfile_path}#{$groupid}_MF_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}#{$groupid}_MF_#{zonalcont}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item == "lat" || item == "plev" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) #{zonalcont}"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys - gphys_con)
      end
    }
  end



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

  def nc_sh_psuv_anm(zonalcont="control")

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

    gphys_ps = ( gphys_ps - gphys_ps_con.mean(0) )
    gphys_ps = gphys_ps.lon_lotate unless $expID =~ /H1998/
    gphys_ps =
      gphys_ps.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")
    gphys_u = @data.gphys_open("ml_u").cut(true,true,100000)
    gphys_v = @data.gphys_open("ml_v").cut(true,true,100000)
    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,true,100000)
    gphys_v_con = @data_con.gphys_open("ml_v").cut(true,true,100000)

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

    unless $expID =~ /H1998/
      gphys_u =  gphys_u.lon_lotate
      gphys_v =  gphys_v.lon_lotate
    end

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

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

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

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

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

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

    gphys_phi = ( gphys_phi - gphys_phi_con.mean(0) )
    gphys_phi = gphys_phi.lon_lotate  unless $expID =~ /H1998/
    gphys_phi = gphys_phi.
      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)")

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

    if $expID =~ /H1998/
      gphys_u = ( gphys_u - gphys_u_con.mean(0) )
      gphys_v = ( gphys_v - gphys_v_con.mean(0) )
    else
      gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
      gphys_v = ( gphys_v - gphys_v_con.mean(0) ).lon_lotate
    end

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

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

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

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

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

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

    gphys_phi = ( gphys_phi - gphys_phi_con.mean(0) )
    gphys_phi = gphys_phi.lon_lotate   unless $expID =~ /H1998/
    gphys_phi = gphys_phi.
      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)")
    
    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)

    if $expID =~ /H1998/      
      gphys_u = ( gphys_u - gphys_u_con.mean(0) )
      gphys_v = ( gphys_v - gphys_v_con.mean(0) )
    else
      gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
      gphys_v = ( gphys_v - gphys_v_con.mean(0) ).lon_lotate
    end

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

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

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

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

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

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

    gphys_t = ( gphys_t - gphys_t_con.mean(0) )
    gphys_t = gphys_t.lon_lotate  unless $expID =~ /H1998/
    gphys_t = gphys_t.
      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)")

    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)

    if $expID =~ /H1998/
      gphys_u = ( gphys_u - gphys_u_con.mean(0) )
      gphys_v = ( gphys_v - gphys_v_con.mean(0) )
    else
      gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
      gphys_v = ( gphys_v - gphys_v_con.mean(0) ).lon_lotate
    end

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

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

    mkfig_plot(gphys_t,u_data,v_data)
    
  end
  
  def nc_ml_tuom_anm(zonalcont="control")

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

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

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

    gphys_t = ( gphys_t - gphys_t_con.mean(0) )
    gphys_t = gphys_t.lon_lotate unless $expID =~ /H1998/
    gphys_t = gphys_t.
      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").cut(true,0,true)
    gphys_om = @data.gphys_open("ml_om").cut(true,0,true)
    gphys_u_con = @data_con.gphys_open("ml_u").cut(true,0,true)
    gphys_om_con = @data_con.gphys_open("ml_om").cut(true,0,true)

    if  $expID =~ /H1998/
      gphys_u = ( gphys_u - gphys_u_con.mean(0) )
      gphys_om = ( gphys_om - gphys_om_con.mean(0) )
    else
      gphys_u = ( gphys_u - gphys_u_con.mean(0) ).lon_lotate
      gphys_om = ( gphys_om - gphys_om_con.mean(0) ).lon_lotate
    end

#    gphys_om = gphys_om/100

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

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



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

  # #{$groupid}_GT_control.nc
  def nc_gt_comparing_expid
#    id = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
    id = ["control","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1","flat3keq","Qobs3keq","H1998con","H1998pa","HS1986"]

    $file_label = "#{$rezol}"
    @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item == "time"  || item == "gt_slwu" # gt_slwu は一定値
	id.each{ |expID|
	  $expID = expID
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_GT_#{$expID}.nc")
	  gphys_ary.push(@data.gphys_open(item).
			 set_att("line_name",$expID) )
	}
	mkfig_plot(gphys_ary)
      end
    }

  end

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal_comparing_expid
#    id = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
    id = ["control","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1","flat3keq","Qobs3keq","H1998con","H1998pa","HS1986"]
    $file_label = "#{$rezol}"
    @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      gphys_ary = []
      unless item == "lon" || item == "lat" 
	id.each{ |expID|
	  $expID = expID
	  @data = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc")
	  gphys = @data.gphys_open(item).mean(0)
  	  gphys_ary.push(gphys.rename("#{gphys.name}_zonal").
			 set_att("line_name",$expID) )
	}
	mkfig_plot(gphys_ary)
      end
    }
  end

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_anm_comparing_expid
    $file_label = "#{$rezol}"

    expID_con_hash =  {
      "1keq" => "control", 
      "3keq" => "control", 
      "3kw1" => "control", 
      "flat3keq" => "flat",
      "Qobs3keq" => "Qobs",
      "H1998pa" => "H1998con"
    }
    id     = ["3keq","flat3keq","Qobs3keq","H1998pa"]

    gphys_ary_tppn = []
    gphys_ary_slh = []
    id.size.times{ |num|
      @data_con = 
	Ape.new("#{$ncfile_path}#{$groupid}_SH_#{expID_con_hash[id[num]]}.nc")
      @data     = Ape.new("#{$ncfile_path}#{$groupid}_SH_#{id[num]}.nc")

      # 赤道上 降水量 (tppn) 
      gphys_con = @data_con.gphys_open("sh_tppn")
      gphys     = @data.gphys_open("sh_tppn")
      if $expID == "H1998con" || $expID == "H1998pa"
	gphys = (gphys - gphys_con.mean(0))
      else
	gphys = (gphys - gphys_con.mean(0)).lon_lotate
      end
      gphys =  gphys.cut(true,0).rename("sh_tppn_lat0_anm"). 
	set_att("ape_name", "total_precipitation_flux"). 
	set_att("line_name","#{expID_con_hash[id[num]]}-#{id[num]}")
      gphys_ary_tppn.push(gphys)
      
      # 赤道上 蒸発量 (slh) 
      gphys = @data.gphys_open("sh_slh")
      gphys_con = @data_con.gphys_open("sh_slh")
      if $expID == "H1998con" || $expID == "H1998pa"
	gphys = (gphys - gphys_con.mean(0))
      else
	gphys = (gphys - gphys_con.mean(0)).lon_lotate
      end
      gphys = gphys.cut(true,0).rename("sh_slh_lat0_anm").
	set_att("line_name","#{expID_con_hash[id[num]]}-#{id[num]}")
      gphys_ary_slh.push(gphys)

    }
    mkfig_plot(gphys_ary_tppn)
    mkfig_plot(gphys_ary_slh)

  end

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_ml_zonal_comparing_expid
#    id = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
    id = ["control","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1","flat3keq","Qobs3keq","H1998con","H1998pa","HS1986"]
    $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}_SH_#{$expID}.nc"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      gphys = @data.gphys_open("ml_u").cut(true,0,true).mean(0)
      gphys_ary.push(gphys.rename("#{gphys.name}_lat0_zonal").
		     set_att("line_name",$expID) )
    }
    mkfig_plot(gphys_ary)

    # MSE (moist static energy) = CpT + gz + Lq
    # DSE (dry static energy) = CpT + gz
    gphys_dse_ary = []
    gphys_mse_ary = []
    id.each{ |expID|
      $expID = expID
      $groupid = $groupid_hash[$rezol]
      $ncfile_path = "/home/yukiko/work/ape/NetCDF/#{$rezol}/"
      puts "#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc"
      @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      gphys_t = @data.gphys_open("ml_t").mean(0).cut(20,true)
      gphys_q = @data.gphys_open("ml_q").mean(0).cut(20,true)
      gphys_phi = @data.gphys_open("ml_phi").mean(0).cut(20,true)*9.8
      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}_GT_control.nc
  def nc_gt_comparing_rezol
    $file_label = "AGCM5_#{$expID}"
    rezol = ["AGCM5_adj", "AGCM5_kuo", "AGCM5_kuo-nosc"]
    $groupid = $groupid_hash[rezol[0]]
    @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
	  $groupid = $groupid_hash[$rezol]
	  $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
    }

  end

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_sh_zonal_comparing_rezol
    $file_label = "AGCM5_#{$expID}"
    rezol = ["AGCM5_adj", "AGCM5_kuo", "AGCM5_kuo-nosc"]
    $groupid = $groupid_hash[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" 
	rezol.each{ |rezol_item|
	  $rezol = rezol_item
	  $groupid = $groupid_hash[$rezol]
	  $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) )
	  puts $rezol
	}
	mkfig_plot(gphys_ary)
      end
    }
  end

  #  #{$groupid}_SH_control.nc@zonal_mean
  def nc_ml_zonal_comparing_rezol
    $file_label = "AGCM5_#{$expID}"
    rezol = ["AGCM5_adj", "AGCM5_kuo", "AGCM5_kuo-nosc"]
    # U の赤道鉛直断面
    gphys_ary = []
    rezol.each{ |rezol_item|
      $rezol = rezol_item
      $ncfile_path = "/home/yukiko/work/ape/NetCDF/#{$rezol}/"
      puts $ncfile_path
      $groupid = $groupid_hash[$rezol]
      @data = Ape.new("#{$ncfile_path}#{$groupid}_ML_#{$expID}.nc")
      gphys = @data.gphys_open("ml_u").cut(true,0,true).mean(0)
      gphys_ary.push(gphys.rename("#{gphys.name}_lat0_zonal").
		     set_att("line_name",$rezol) )
    }
    mkfig_plot(gphys_ary)

    # 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]
      @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)*9.8
       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

    $file_label = "AGCM5_#{$expID}"
    rezol = ["AGCM5_adj", "AGCM5_kuo", "AGCM5_kuo-nosc"]

    expID_con_hash =  {
      "1keq" => "control", 
      "3keq" => "control", 
      "3kw1" => "control", 
      "flat3keq" => "flat",
      "Qobs3keq" => "Qobs",
      "H1998pa" => "H1998con"
    }
    $file_label = "AGCM5_#{$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
      $groupid = $groupid_hash[$rezol]
      $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) 
      gphys_con = @data_con.gphys_open("sh_tppn")
      gphys     = @data.gphys_open("sh_tppn")
      if $expID == "H1998con" || $expID == "H1998pa"
	gphys = (gphys - gphys_con.mean(0))
      else
	gphys = (gphys - gphys_con.mean(0)).lon_lotate
      end
      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")
      if $expID == "H1998con" || $expID == "H1998pa"
	gphys = (gphys - gphys_con.mean(0))
      else
	gphys = (gphys - gphys_con.mean(0)).lon_lotate
      end
      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

end