#!/usr/bin/env ruby
# ----------------------------------------------
# local load path

# $local_path = '/work11/ape/yukiko/lib'
# $local_path = '/home/yukiko/eva01/work11/ape/yukiko/lib'
 $local_path = '/home/yukiko/work/ape/yukiko/lib'
# $local_path = '/home/yukiko/tmp/ape-data/lib'
$: << $local_path

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

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

END{

host = "eva01"
 a = Ape_mkfig.new 3
# a = Ape_mkfig.new 2
# a = Ape_mkfig.new 1

#  ["con","mradlAa","mradlAb","mradlAc","mradlAd"].each{ |item|
  ["con","a","b","c","d"].each{ |item|
#    $expid = item
#    $cum = "adj"
#    set_path
#    a.nc_tr_agcm5
#    a.nc_tr_wkplot_agcm5

    $expid = item
    $cum = "kuo"
    set_path
    a.nc_tr_agcm5
#    a.nc_tr_wkplot_agcm5
  }


}

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


def set_path
#  $fig_path = "/work11/ape/yukiko/figs/tmp/"
#  $fig_path = "/home/yukiko/eva01/home/yukiko/work/aqua-mradl-ronbun/figs/tmp/"
  $fig_path = "/home/yukiko/work/aqua3d/yukiko/sh/ruby/mradlstdexp/mlshtr/figs/"
#  $ncfile_path = "/work11/aqua3d/yukiko/mradl/data/wk-plot-data/"
  $ncfile_path = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/"
#  $file_label = "AGCM5_#{$cum}_#{$expid}"

  if $expid == "con"
    $file_label = "#{$cum}_control"
#    $file_label = "#{$cum}_con"
  else
    $file_label = "#{$cum}_#{$expid}"
    $expid = "mradlA" + $expid
  end

  $groupid = "AGCM5-#{$cum}-#{$expid}"
end

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

class Ape_mkfig
  # #{$groupid}_TR_control.nc
  def nc_tr_agcm5

    netcdf_var = 
      { "rain_sum"   => "rain",
      "u_out_z12" => "u",
      "u_out_z13" => "u",
      "v_out_z12" => "v", 
      "v_out_z13" => "v",
      "gph_out_z12" => "z",
      "gph_out_z13" => "z", 
     "t_out_z13" => "t",
     "t_out_z12" => "t",
     "t_out_z8" => "t",
     "q_out_z6" => "q", 
     "q_out_z5" => "q", 
     "q_out_z3" => "q", 
      "ps_out"    => "ps"  }

    netcdf_var_rename = 
      { "rain_sum"   => "tr_tppn",
      "u_out_z12" => "tr_u_z12",
      "u_out_z13" => "tr_u_z13",
      "v_out_z12" => "tr_v_z12", 
      "v_out_z13" => "tr_v_z13",
      "gph_out_z12" => "tr_gph_z12",
      "gph_out_z13" => "tr_gph_z13", 
     "t_out_z13"  => "t_out_z13",
     "t_out_z12"  => "t_out_z12",
     "t_out_z8"   => "t_out_z8",
     "q_out_z6"   => "q_out_z6", 
     "q_out_z5"   => "q_out_z5", 
     "q_out_z3"   => "q_out_z3", 
      "ps_out"    => "tr_mslp"  }

    netcdf_apename = 
      { "rain_sum"   => "precipitation_flux",
      "u_out_z12" => "eastward_wind",
      "u_out_z13" => "eastward_wind",
      "v_out_z12" => "northward_wind",
      "v_out_z13" => "northward_wind",
      "gph_out_z12" => "geopotential_height",
      "gph_out_z13" => "geopotential_height", 
      "t_out_z13"  => "temperature",
      "t_out_z12"  => "temperature",
      "t_out_z8"   => "temperature",
      "q_out_z6"   => "specific_humidity", 
      "q_out_z5"   => "specific_humidity", 
      "q_out_z3"   => "specific_humidity", 
      "ps_out"    => "air_pressure_at_sea_level"  }

    netcdf_units = 
      { "rain_sum"   => "kg m-2 s-1",
      "u_out_z12" => "m s-1",
      "u_out_z13" => "m s-1",
      "v_out_z12" => "m s-1", 
      "v_out_z13" => "m s-1",
      "gph_out_z12" => "m",
      "gph_out_z13" => "m", 
      "t_out_z13"  => "K",
      "t_out_z12"  => "K",
      "t_out_z8"   => "K",
      "q_out_z6"   => "1", 
      "q_out_z5"   => "1", 
      "q_out_z3"   => "1", 
      "ps_out"    => "Pa"  }

    lat_lost_axis = ["latitude=1.395307"]

#    ["rain_sum","u_out_z12","u_out_z13","v_out_z12","v_out_z13","gph_out_z12","gph_out_z13","ps_out"].each { |varname|
    ["rain_sum"].each { |varname|

      @data = Ape.new("#{$ncfile_path}/#{$cum}700#{$expid}400/#{varname}.nc")

      name = netcdf_var_rename[varname]
      
      if name =~ /mslp|z13|z12/
	if name =~ /mslp/
	  lost_axis = []
	elsif name =~ /z13/
	  lost_axis = ["sigma=0.17457"]
	elsif name =~ /z12/
	  lost_axis = ["sigma=0.22953"]
	end

	# x-t ダイヤグラム
	gphys = @data.gphys_open(netcdf_var[varname]).rename(name)
	gphys =  gphys.cut(true,0,true)[true,-400..-1]

	gphys = gphys.
	  set_att("ape_name",netcdf_apename[varname]).
	  set_att("units",netcdf_units[varname])

	if varname == "ps_out"
	  gphys = gphys * 100.0
	end

	grid_1 =
	  Axis.new().
	  set_pos(VArray.new(NArray.sfloat(400).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 )

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

#	mkfig_plot((gphys - gphys.mean(0)).set_lost_axes(""))

	# 時空間スペクトル
	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(lat_lost_axis).
	  add_lost_axes(lost_axis)
	
	dim1 = gphys.coord(0).shape.to_s.to_i
	dim2 = gphys.coord(1).shape.to_s.to_i
	mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..80] )
	
      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open(netcdf_var[varname]).rename(name)
#	gphys = @data.gphys_open("tr_tppn")
	gphys = gphys.cut(true,0,true)[true,-400..-1]*0.4*10**(-6)

	gphys = gphys.
	  set_att("ape_name",netcdf_apename[varname]).
	  set_att("units",netcdf_units[varname])

	grid_1 =
	  Axis.new().
	  set_pos(VArray.new(NArray.sfloat(400).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 )

	mkfig_plot(gphys.set_lost_axes(lat_lost_axis))
#	mkfig_plot(gphys.rename("tr_tppn_mono").set_lost_axes(""))
	
      end
    }
  end

end


class Ape

  
  def mkdump(gphys, gphys_u=nil, gphys_v=nil)
    gropn(3) if $gropn == nil

    if gphys_u == nil || gphys_v == nil then
      plot(gphys)
    else
      plot(gphys, gphys_u, gphys_v)
    end
    
    gphys = gphys[0]     if gphys.class == Array
    
    file_name = "#{$file_label.gsub("control","con")}_#{gphys.name}.gif"
    `rm dcl_* tmp.pnm`  
    plot(gphys)
    `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm;done`
#    `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 840 630 > tmp.pnm;done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm`  

    $file_name = $pre_file

  end


  # トーン, コンターパターンの設定 (変数固有)
  def plot_set(gphys)

    # コンターを描かない条件
    if gphys.name =~ /tr_/ || gphys.name =~ /horinout/  
      $cont_flag = false 
    else
      $cont_flag = true
    end

    # トーンを描かない条件
    if gphys.name == "comp_point"
      $tone_flag = false 
    else
      $tone_flag = true
    end

    # 初期化
    $fig_set_hash = {"window" => nil}
    $line_hash = { 'exchange'=>false, 'index'=>3 }


    if gphys.name == "tr_tppn" 
      # 塚原カラー

#      patterns = NArray[35999, 45999, 55999, 70999, 75999, 85999]
#      levels = NArray[10.0**(-6), 10.0**(-5.5), 10.0**(-5), 10.0**(-4.5), 10.0**(-4), 10.0**(-3.5), DCL.glpget('rmiss')]

      patterns = NArray[ 45999, 75999]
      levels = NArray[(2.3)*10**(-5), (5.8)*10**(-5), DCL.glpget('rmiss')]

      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'int'=>0.5}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1
      }

    elsif gphys.name == "tr_u_z13" || gphys.name == "tr_u_z12" || gphys.name == "tr_v_z13" || gphys.name == "tr_v_z12" 
      levels = NArray[DCL.glpget('rmiss'), -10, -5, 0, 5, 10, 20, DCL.glpget('rmiss')]
      patterns = NArray[30999, 35999, 40999, 55999, 70999, 75999, 85999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>0,
	"tick1"=>20,"tick2"=>1
      }

    elsif gphys.name == "tr_gph_z13" || gphys.name == "tr_gph_z12" 
      levels = NArray[DCL.glpget('rmiss'), -30, -10, 0, 10, 30, 50, DCL.glpget('rmiss')]
      patterns = NArray[30999, 35999, 40999, 55999, 70999, 75999, 85999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>0,
	"tick1"=>20,"tick2"=>1
      }

    elsif gphys.name == "tr_mslp" 
      levels = NArray[DCL.glpget('rmiss'), -450, -300, -150, 0, 150, 300, DCL.glpget('rmiss')]
      patterns = NArray[30999, 35999, 40999, 55999, 70999, 75999, 85999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>0,
	"tick1"=>20,"tick2"=>1
      }

    elsif gphys.name == "tr_tppn_spct" 
      levels = NArray[1e-12, 1e-11, 1e-10, 1e-9, 1e-8 ]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'int'=>1e-10}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    elsif gphys.name == "tr_u_z13_spct" || gphys.name == "tr_u_z12_spct" 
      levels = NArray[1e-3, 1e-2, 1e-1, 1, 10 ]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'int'=>1e-1}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    elsif gphys.name == "tr_v_z13_spct" || gphys.name == "tr_v_z12_spct" 
      levels = NArray[1e-3, 1e-2, 1e-1, 1, 10 ]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'int'=>1e-1}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1, 
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    elsif gphys.name == "tr_gph_z13_spct" || gphys.name == "tr_gph_z12_spct" 
      levels = NArray[1e-3, 1e-2, 1e-1, 1, 10 ]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'int'=>1e-1}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1, 
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    elsif gphys.name == "tr_mslp_spct" 
      levels = NArray[1e+0, 1e+1, 1e+2, 1e+3, 1e+4 ]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'int'=>1e+2}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

      # W-K plot
    elsif gphys.name == "tr_tppn_sym_spct" || gphys.name == "tr_u_z13_sym_spct" || gphys.name == "tr_u_z12_sym_spct" ||  gphys.name == "tr_v_z13_sym_spct" ||  gphys.name == "tr_v_z12_sym_spct" || gphys.name == "tr_gph_z13_sym_spct" || gphys.name == "tr_gph_z12_sym_spct"|| gphys.name == "tr_mslp_sym_spct" 
      levels = NArray[1.2, 1.4, 1.6, 2.0, DCL.glpget('rmiss')]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    elsif gphys.name == "tr_tppn_asym_spct" || gphys.name == "tr_u_z13_asym_spct" || gphys.name == "tr_u_z12_asym_spct" ||  gphys.name == "tr_v_z13_asym_spct" ||  gphys.name == "tr_v_z12_asym_spct"||  gphys.name == "tr_gph_z13_asym_spct" ||  gphys.name == "tr_gph_z12_asym_spct" || gphys.name == "tr_mslp_asym_spct" 
      levels = NArray[0.1, 0.3, 0.5, 0.7, DCL.glpget('rmiss')]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    else
      plot_def(gphys)
    end

  end

end