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

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

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

load "#{$local_path}/ape-view.rb"
load "/home/yukiko/aqua-mradl-ronbun/tools/work-AGCM5-tr-tq-plot.rb"

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

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

# ----------------------------------------------
groupid = ["aguforape","agcm5_adj","agcm5_kuo"]
host = "eva01"
set_dir_id  = groupid[0]

END{

  ["b","a","c","d","con"].each{ |item|
    $expid = item
    $cum = "adj"
    set_path
    a.nc_tr_spctfilt
    a.nc_tr_spctfilt_overplot
    $expid = item
    $cum = "kuo"
    set_path
    a.nc_tr_spctfilt
    a.nc_tr_spctfilt_overplot    
  }

}

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


def set_path
  $fig_path = "/home/yukiko/aqua-mradl-ronbun/figs/tmp/"
  $ncfile_path = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/spctfilt/"
  #  $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

  def nc_tr_spctfilt
    
    filt_array = ["ad","k","g","n"]
    filt_hash =  { 
      "ad" => "advect", 
      "k"  => "kelvin",
      "g"  => "grav", 
      "n"  => "no", 
    }

    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}_#{filt_hash[filter]}filt.nc")

      p "#{dir_path}AGCM5-#{$cum}-#{$expid}_#{filt_hash[filter]}filt.nc"

      # PR 
      if filter == "n"
        gphys = @data.gphys_open("tr_tppn")
      else
        gphys = @data.gphys_open("tr_tppn_#{filter}filt")
      end
      gphys = gphys.add_lost_axes( gphys.get_att("lost_axis") )
      gphys = gphys[true,-401..-1]
      mkfig_plot(gphys)
      spct_plot(gphys)
      
      # T (0.57->0.54945)
      if filter == "n"
        gphys = @data.gphys_open("tr_t").
          rename("tr_t_z8")
      else
        gphys = @data.gphys_open("tr_t_#{filter}filt"). 
          rename("tr_t_z8_#{filter}filt")
      end
      gphys = gphys.cut(true,0.54945,true)
      gphys = gphys.
        set_lost_axes(gphys.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","")).
        add_lost_axes(gphys.get_att("lost_axis") )
      gphys = gphys[true,-401..-1]
      mkfig_plot((gphys - gphys.mean(0)).
                   add_lost_axes("(diff) from (mean) zonal"))
      spct_plot(gphys)
      
      # Q (0.85->0.82977)
      if filter == "n"
        gphys = @data.gphys_open("tr_q").
          rename("tr_q_z5")
      else
        gphys = @data.gphys_open("tr_q_#{filter}filt"). 
          rename("tr_q_z5_#{filter}filt")
      end
      gphys = gphys.cut(true,0.82977,true)
      gphys = gphys.
        set_lost_axes(gphys.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","")).
        add_lost_axes(gphys.get_att("lost_axis") )
      gphys = gphys[true,-401..-1]
      mkfig_plot((gphys - gphys.mean(0)).
                   add_lost_axes("(diff) from (mean) zonal"))
      spct_plot(gphys)
      
      # U (0.25->0.22953)
      if filter == "n"
        gphys = @data.gphys_open("tr_u").
          rename("tr_u_z12")
      else
        gphys = @data.gphys_open("tr_u_#{filter}filt"). 
          rename("tr_u_z12_#{filter}filt")
      end
      gphys = gphys.cut(true,0.22953,true)
      gphys = gphys.
        set_lost_axes(gphys.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","")).
        add_lost_axes(gphys.get_att("lost_axis") )
      gphys = gphys[true,-401..-1]
      mkfig_plot((gphys - gphys.mean(0)).
                   add_lost_axes("(diff) from (mean) zonal"))
      spct_plot(gphys)
    }    

  end
    
  def spct_plot(gphys)
    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")})2").
      add_lost_axes(gphys.lost_axes)
    dim1 = gphys.coord(0).shape.to_s.to_i
    dim2 = gphys.coord(1).shape.to_s.to_i
    mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
  end

  def nc_tr_spctfilt_overplot
    
    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}_#{filt_hash[filter]}filt.nc")
      @data_filt =
        Ape.new("#{dir_path}AGCM5-#{$cum}-#{$expid}_nofilt.nc")
      
      # PR 
      gphys = @data_filt.gphys_open("tr_tppn")
      filt = @data.gphys_open("tr_tppn_#{filter}filt")
      filt = filt.rename("tr_tppn_#{filter}filt_op")
      gphys = gphys.rename("tr_tppn_#{filter}filt_op")
      gphys = gphys.add_lost_axes( gphys.get_att("lost_axis") )
      mkfig_plot_comp(gphys, filt)
      
      # T (0.57->0.54945)
      gphys =
        @data_filt.gphys_open("tr_t").cut(true,0.54945,true)
      filt = 
        @data.gphys_open("tr_t_#{filter}filt").cut(true,0.54945,true)
      gphys = gphys.rename("tr_t_z8_#{filter}filt_op")
      filt = filt.rename("tr_t_z8_#{filter}filt_op")
      gphys = gphys.
        set_lost_axes(gphys.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","")).
        add_lost_axes(gphys.get_att("lost_axis") )
      gphys = (gphys - gphys.mean(0)).
        add_lost_axes("(diff) from (mean) zonal")
      filt = (filt - filt.mean(0))
      mkfig_plot_comp(gphys, filt)

      # Q (0.85->0.82977)
      gphys =
        @data_filt.gphys_open("tr_q").cut(true,0.82977,true)
      filt = 
        @data.gphys_open("tr_q_#{filter}filt").cut(true,0.82977,true)
      gphys = gphys.rename("tr_q_z5_#{filter}filt_op")
      filt = filt.rename("tr_q_z5_#{filter}filt_op")
      gphys = gphys.
        set_lost_axes(gphys.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","")).
        add_lost_axes(gphys.get_att("lost_axis") )
      gphys = (gphys - gphys.mean(0)).
        add_lost_axes("(diff) from (mean) zonal")
      filt = (filt - filt.mean(0))
      mkfig_plot_comp(gphys, filt)
      
      # U (0.25->0.22953)
      gphys =
        @data_filt.gphys_open("tr_u").cut(true,0.22953,true)
      filt = 
        @data.gphys_open("tr_u_#{filter}filt").cut(true,0.22953,true)
      gphys = gphys.rename("tr_u_z12_#{filter}filt_op")
      filt = filt.rename("tr_u_z12_#{filter}filt_op")
      gphys = gphys.
        set_lost_axes(gphys.lost_axes.to_s.gsub("sigmap","sigma").gsub(" 1","")).
        add_lost_axes(gphys.get_att("lost_axis") )
      gphys = (gphys - gphys.mean(0)).
        add_lost_axes("(diff) from (mean) zonal")
      filt = (filt - filt.mean(0))
      mkfig_plot_comp(gphys, filt)
      
    }    

  end


end