#!/usr/bin/env ruby

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

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

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

load "#{$local_path}/ape-view.rb"
load "#{$local_path}/lib-ape-view.rb"
load "#{$local_path}/lib-ape-base.rb"
load "/home/yukiko/work/aqua3d/yukiko/sh/ruby/mradlstdexp/ronbun/work-AGCM5-tr-tq-plot.rb"

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


END{

  a = Ape_mkfig.new 3

#  ["c", "con", "a", "b", "d" ].each{ |expid|
  ["c" ].each{ |expid|
    $expid = expid
    $cum = "kuo"
    set_path
#    a.gr_comp_tuv_horizontal
#    a.gr_comp_tuv
#    a.gr_comp_point
    $expid = expid
    $cum = "adj"
    set_path
    a.gr_comp_tuv_horizontal
#    a.gr_comp_tuv
#    a.gr_comp_point

  }
}

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


def set_path
  $fig_path = "/home/yukiko/work/aqua3d/yukiko/sh/ruby/mradlstdexp/ronbun/figs/"
  $ncfile_path = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/util/work-AGCM5-spctfilt/spctfilt/"
  if $expid == "con"
    $file_label = "#{$cum}_control"
  else
    $file_label = "#{$cum}_#{$expid}"
    $expid = "mradlA" + $expid
  end
  $groupid = "AGCM5-#{$cum}-#{$expid}"
end

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

class Ape_mkfig

  def mkfig_plot_comp(gphys, gphys_comp)

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

  end

  def gr_comp_point(compmethod=nil)
    
    if compmethod == nil 
      compmethod="threshold"
    end
    
    filt_array = ["ad","k","g"]
    filt_hash =  { 
      "ad" => "advect", 
      "k"  => "kelvin",
      "g"  => "grav"
    }

    filt_array.each{ |filter|

      dir_path = 
        "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/spctfilt/"
      @data =
        Ape.new("#{dir_path}AGCM5-#{$cum}-#{$expid}_nofilt.nc")
      
      # PR 
      rain = @data.gphys_open("tr_tppn")
      rain = rain.add_lost_axes(rain.get_att("lost_axis"))
      rain = rain[true,-401..-1]
      rain = rain.
        rename("comp_point_tr_tppn_#{filter}filtreal")

      # U (0.25->0.22953)
      u250 = @data.gphys_open("tr_u").cut(true,0.22953,true)
      u250= u250.
        set_lost_axes(u250.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","") ).
        add_lost_axes(u250.get_att("lost_axis"))
      u250 = u250[true,-401..-1]
      u250 = ( u250 - u250.mean(0) )
      u250 = u250.
        add_lost_axes("(diff) from (mean) zonal"). 
        rename("comp_point_tr_u_z12_#{filter}filtreal")
      
      # T (0.57->0.54945)
      t570 = @data.gphys_open("tr_t").cut(true,0.54945,true)
      t570= t570.
        set_lost_axes(t570.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","") ).
        add_lost_axes(t570.get_att("lost_axis"))
      t570 = t570[true,-401..-1]
      t570 = ( t570 - t570.mean(0) )
      t570 = t570.
        add_lost_axes("(diff) from (mean) zonal"). 
        rename("comp_point_tr_t_z8_#{filter}filtreal")
      
      # Q (0.85->0.82977)
      q850 = @data.gphys_open("tr_q").cut(true,0.82977,true)
      q850= q850.
        set_lost_axes(q850.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","")).
        add_lost_axes(q850.get_att("lost_axis"))
      q850 = q850[true,-401..-1]
      q850 = ( q850 - q850.mean(0) )
      q850 = q850.
        add_lost_axes("(diff) from (mean) zonal"). 
        rename("comp_point_tr_q_z5_#{filter}filtreal")

      comp_dir_path = 
        "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/util/work-AGCM5-spctfilt/spctfilt/"        
      file_name = 
        "#{comp_dir_path}AGCM5-#{$cum}-#{$expid}_#{filt_hash[filter]}filtreal-composite.nc"

      @data = Ape.new(file_name)
      comp = @data.gphys_open("comp_point")
      
      lost_axis = 
        comp.data.get_att("lost_axis").to_s.split("\n")[-2..-1]
      
      rain = rain.add_lost_axes(lost_axis )
      u250 = u250.add_lost_axes(lost_axis )
      q850 = q850.add_lost_axes(lost_axis )
      t570 = t570.add_lost_axes(lost_axis )

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

  end
  
  def gr_comp_tuv(compmethod=nil)
    
    if compmethod == nil 
      compmethod="threshold"
    end
    
    filt_array = ["ad","k","g"]
    filt_hash =  { 
      "ad" => "advect", 
      "k"  => "kelvin",
      "g"  => "grav"
    }

    filt_array.each{ |filter|
      
      comp_dir_path = 
        "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/util/work-AGCM5-spctfilt/spctfilt/"        
      file_name = 
        "#{comp_dir_path}AGCM5-#{$cum}-#{$expid}_#{filt_hash[filter]}filtreal-composite.nc"
    
      @data = Ape.new(file_name)
      t_comp = @data.gphys_open("comp_t_#{filter}filtreal").
        rename("comp_tuom_#{filter}filtreal")
      q_comp = @data.gphys_open("comp_q_#{filter}filtreal").
        rename("comp_quom_#{filter}filtreal")
      tconv_comp = @data.gphys_open("comp_tconv_#{filter}filtreal").
        rename("comp_tconvuom_#{filter}filtreal") 
      u_comp = @data.gphys_open("comp_u_#{filter}filtreal").
        rename("comp_uuom_#{filter}filtreal") 
      om_comp = @data.gphys_open("comp_om_#{filter}filtreal")[true,1..-1]

      lost_axes_diff = t_comp.get_att("lost_axis").to_s.split("\n")
      lost_axes = tconv_comp.get_att("lost_axis").to_s.split("\n")
      
      t_comp = t_comp.set_lost_axes(lost_axes_diff)
      u_comp = u_comp.set_lost_axes(lost_axes_diff)
      q_comp = q_comp.set_lost_axes(lost_axes_diff)
      tconv_comp = tconv_comp.set_lost_axes(lost_axes)
      
      # ベクトルの間引
      u_data  = NArray.sfloat( u_comp.axis(0).pos.val.size , 
                               u_comp.axis(1).pos.val.size ).fill!(0.0)
      
      om_data = NArray.sfloat( om_comp.axis(0).pos.val.size , 
                               om_comp.axis(1).pos.val.size ).fill!(0.0)
      
      if compmethod =~ /_large/ 
        num_l = u_comp.axis(0).pos.val.size/60
      else
        num_l = u_comp.axis(0).pos.val.size/120
      end
      num_z = u_comp.axis(1).pos.val.size/48
    
      u_comp.axis(0).pos.val.size.times{ |num|
        if (num % num_l) == 0
	  u_data[num,true]  =  u_comp.data.val[num,true]
  	  om_data[num,true] =  om_comp.data.val[num,true]
        elsif
	  u_data[num,true]  = 0.0
	  om_data[num,true] = 0.0
        end
      }
    
      num_z =1 if num_z < 1

      u_comp.axis(1).pos.val.size.times{ |num|
        if (num % num_z) == 0
	  u_data[true,num]  =  u_data[true,num]
  	  om_data[true,num] =  om_data[true,num]
        elsif
	  u_data[true,num]  = 0.0
	  om_data[true,num] = 0.0
        end
      }

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

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

        mkfig_plot(u_comp,
	           u_data,
	           -om_data )

        mkfig_plot(q_comp, 
	           u_data, 
	           -om_data )

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

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

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

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

     }

  end
  
  def gr_comp_tuv_horizontal(compmethod=nil)
    
    if compmethod == nil 
      compmethod="threshold"
    end
    
    filt_array = ["ad","k","g"]
    filt_hash =  { 
      "ad" => "advect", 
      "k"  => "kelvin",
      "g"  => "grav"
    }

    filt_array.each{ |filter|


      comp_dir_path = 
        "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/util/work-AGCM5-spctfilt/spctfilt/"        
      file_name = 
        "#{comp_dir_path}AGCM5-#{$cum}-#{$expid}_#{filt_hash[filter]}filtreal-composite-horizontal.nc"
    
      @data = Ape.new(file_name)
      tppn_comp = @data.gphys_open("comp_tppn_#{filter}filtreal").
        rename("comp_tppn_#{filter}filtreal")
      psuv_comp = @data.gphys_open("comp_mslp_#{filter}filtreal").
        rename("comp_psuv_#{filter}filtreal")
      phiuv850_comp = @data.gphys_open("comp_gph_z5_#{filter}filtreal").
        rename("comp_gphuv_z5_#{filter}filtreal")
      phiuv250_comp = @data.gphys_open("comp_gph_z12_#{filter}filtreal").
        rename("comp_gphuv_z12_#{filter}filtreal")
      tuv570_comp = @data.gphys_open("comp_t_z8_#{filter}filtreal").
        rename("comp_tuv_z8_#{filter}filtreal")
      quv850_comp = @data.gphys_open("comp_q_z5_#{filter}filtreal").
        rename("comp_quv_z5_#{filter}filtreal")
      
      u250s_comp = @data.gphys_open("comp_u_z12_#{filter}filtreal")
      u570s_comp = @data.gphys_open("comp_u_z8_#{filter}filtreal")
      u850s_comp = @data.gphys_open("comp_u_z5_#{filter}filtreal")
      u950s_comp = @data.gphys_open("comp_u_z3_#{filter}filtreal")
      v250s_comp = @data.gphys_open("comp_v_z12_#{filter}filtreal")
      v570s_comp = @data.gphys_open("comp_v_z8_#{filter}filtreal")
      v850s_comp = @data.gphys_open("comp_v_z5_#{filter}filtreal")
      v950s_comp = @data.gphys_open("comp_v_z3_#{filter}filtreal")
      
      tppn_comp = tppn_comp.
        set_lost_axes( tppn_comp.get_att("lost_axis").split(/\n/) )
      psuv_comp = psuv_comp.
        set_lost_axes( psuv_comp.get_att("lost_axis").split(/\n/) )
      phiuv850_comp = phiuv850_comp.
        set_lost_axes( phiuv850_comp.get_att("lost_axis").split(/\n/) )
      phiuv250_comp = phiuv250_comp.
        set_lost_axes( phiuv250_comp.get_att("lost_axis").split(/\n/) )
      tuv570_comp = tuv570_comp.
        set_lost_axes( tuv570_comp.get_att("lost_axis").split(/\n/) )
      quv850_comp = quv850_comp.
        set_lost_axes( quv850_comp.get_att("lost_axis").split(/\n/) )
      
      # ベクトルの間引
      u250s_data = NArray.sfloat( u250s_comp.axis(0).pos.val.size , 
                                  u250s_comp.axis(1).pos.val.size ).fill!(0.0)
      v250s_data = NArray.sfloat( v250s_comp.axis(0).pos.val.size , 
                                  v250s_comp.axis(1).pos.val.size ).fill!(0.0)
      
      u850s_data = NArray.sfloat( u850s_comp.axis(0).pos.val.size , 
                                  u850s_comp.axis(1).pos.val.size ).fill!(0.0)
      v850s_data = NArray.sfloat( v850s_comp.axis(0).pos.val.size , 
                                  v850s_comp.axis(1).pos.val.size ).fill!(0.0)
      
      u950s_data = NArray.sfloat( u950s_comp.axis(0).pos.val.size , 
                                  u950s_comp.axis(1).pos.val.size ).fill!(0.0)
      v950s_data = NArray.sfloat( v950s_comp.axis(0).pos.val.size , 
                                  v950s_comp.axis(1).pos.val.size ).fill!(0.0)
      
      u570s_data = NArray.sfloat( u570s_comp.axis(0).pos.val.size , 
                                  u570s_comp.axis(1).pos.val.size ).fill!(0.0)
      v570s_data = NArray.sfloat( v570s_comp.axis(0).pos.val.size , 
                                  v570s_comp.axis(1).pos.val.size ).fill!(0.0)
      
      num_l = u250s_comp.axis(0).pos.val.size/120
      num_m = u250s_comp.axis(1).pos.val.size/60*3
      
      u250s_comp.axis(0).pos.val.size.times{ |num|
        if (num % num_l) == 0
          u250s_data[num,true]  =  u250s_comp.data.val[num,true]
          v250s_data[num,true]  =  v250s_comp.data.val[num,true]
          u850s_data[num,true]  =  u850s_comp.data.val[num,true]
          v850s_data[num,true]  =  v850s_comp.data.val[num,true]
          u950s_data[num,true]  =  u950s_comp.data.val[num,true]
          v950s_data[num,true]  =  v950s_comp.data.val[num,true]
          u570s_data[num,true]  =  u570s_comp.data.val[num,true]
          v570s_data[num,true]  =  v570s_comp.data.val[num,true]
          
        elsif
          u250s_data[num,true]  = 0.0
          v250s_data[num,true]  = 0.0
          u850s_data[num,true]  = 0.0
          v850s_data[num,true]  = 0.0
          u950s_data[num,true]  = 0.0
          v950s_data[num,true]  = 0.0
          u570s_data[num,true]  = 0.0
          v570s_data[num,true]  = 0.0
        end
      }
    
      num_m =1 if num_m < 1
    
      u250s_comp.axis(1).pos.val.size.times{ |num|
        if (num % num_m) == 0
          u250s_data[true,num]  =  u250s_data[true,num]
          v250s_data[true,num]  =  v250s_data[true,num]
          u850s_data[true,num]  =  u850s_data[true,num]
          v850s_data[true,num]  =  v850s_data[true,num]
          u950s_data[true,num]  =  u950s_data[true,num]
          v950s_data[true,num]  =  v950s_data[true,num]
          u570s_data[true,num]  =  u570s_data[true,num]
          v570s_data[true,num]  =  v570s_data[true,num]
        elsif
          u250s_data[true,num]  = 0.0
          v250s_data[true,num]  = 0.0
          u850s_data[true,num]  = 0.0
          v850s_data[true,num]  = 0.0
          u950s_data[true,num]  = 0.0
          v950s_data[true,num]  = 0.0
          u570s_data[true,num]  = 0.0
          v570s_data[true,num]  = 0.0
        end
      }
      
      # 描画
      dim = tppn_comp.axis(0).pos.val.size
      
      mkfig_plot(tppn_comp[(dim/3)..(dim*2/3),true] )
      
      mkfig_plot(psuv_comp[(dim/3)..(dim*2/3),true],
                 u950s_data[(dim/3)..(dim*2/3),true], 
                 v950s_data[(dim/3)..(dim*2/3),true])
      
      mkfig_plot(phiuv850_comp[(dim/3)..(dim*2/3),true],
                 u850s_data[(dim/3)..(dim*2/3),true], 
                 v850s_data[(dim/3)..(dim*2/3),true])
      
      mkfig_plot(phiuv250_comp[(dim/3)..(dim*2/3),true],
                 u250s_data[(dim/3)..(dim*2/3),true], 
                 v250s_data[(dim/3)..(dim*2/3),true])
      
      mkfig_plot(tuv570_comp[(dim/3)..(dim*2/3),true],
                 u570s_data[(dim/3)..(dim*2/3),true], 
                 v570s_data[(dim/3)..(dim*2/3),true])
      
      mkfig_plot(quv850_comp[(dim/3)..(dim*2/3),true],
                 u850s_data[(dim/3)..(dim*2/3),true], 
                 v850s_data[(dim/3)..(dim*2/3),true])
      
    }
    
  end
end



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

    plot_main_comp(gphys, gphys_u)

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

    plot_comp(gphys, gphys_u)

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

    $file_name = $pre_file

  end

  def plot_main_comp(gphys, gphys_u=nil)

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

    GGraph.set_fig($fig_set_hash)

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

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

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

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

  end
end

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