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

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

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

load "#{$local_path}/ape-view.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{

sstid = ["control"]

rezolid = [
  "T39L48_non", "T159L48_non", "T39L48_eml", "T159L48_eml"
# "T39L24_non","T39L96_non", "T79L48_non"
#  "T39L24_eml","T39L96_eml", "T79L48_eml", 
#  "T39L48_ias", "T39L48_ksc", "T39L48_kuo", "T39L48_mca" 
]

  rezolid.each { |item|
    $expID = "control"
    $rezol = item
    set_dir(set_dir_id,host)
#    a.nc_tr_spctfilt
    a.nc_tr_spctfilt_overplot

  }

}

class Ape_mkfig

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

    filt_array.each{ |filter|

      dir_path = "/home/yukiko/work/ape/yukiko/data/spctfilt/"
      @data =
        Ape.new("#{dir_path}#{$rezol}_#{$expID}_#{filt_hash[filter]}filt.nc")

      # PR 
      gphys = @data.gphys_open("tr_tppn_#{filter}filt")
      gphys = gphys.add_lost_axes(gphys.get_att("lost_axis"))
      if $rezol == "T159L48_eml"
        gphys = gphys[true,-801..-401]
      else
        gphys = gphys[true,-401..-1]
      end
      mkfig_plot(gphys.lon_lotate)
      spct_plot(gphys)
      
      # T (0.57)
      gphys = @data.gphys_open("tr_t_#{filter}filt").cut(true,0.57,true)
      gphys = gphys.
        add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma")).
        rename("tr_t570s_#{filter}filt")
      if $rezol == "T159L48_eml"
        gphys = gphys[true,-801..-401]
      else
        gphys = gphys[true,-401..-1]
      end
      mkfig_plot((gphys - gphys.mean(0)).
                   add_lost_axes("(diff) from (mean) zonal").lon_lotate)
      spct_plot(gphys)
      
      # Q (0.85)
      gphys = @data.gphys_open("tr_q_#{filter}filt").cut(true,0.85,true)
      gphys = gphys.
        add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma") ).
        rename("tr_q850s_#{filter}filt")
      if $rezol == "T159L48_eml"
        gphys = gphys[true,-801..-401]
      else
        gphys = gphys[true,-401..-1]
      end
      mkfig_plot((gphys - gphys.mean(0)).
                   add_lost_axes("(diff) from (mean) zonal").lon_lotate)
      spct_plot(gphys)
      
      # U (0.25)
      gphys = @data.gphys_open("tr_u_#{filter}filt").cut(true,0.25,true)
      gphys = gphys.
        add_lost_axes(gphys.get_att("lost_axis").gsub("z","sigma") ).
        rename("tr_u250s_#{filter}filt")
      if $rezol == "T159L48_eml"
        gphys = gphys[true,-801..-401]
      else
        gphys = gphys[true,-401..-1]
      end
      mkfig_plot((gphys - gphys.mean(0)).
                   add_lost_axes("(diff) from (mean) zonal").lon_lotate)
      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"
    }

#    DCL.sgscmn(4) 

    filt_array.each{ |filter|

      dir_path = "/home/yukiko/work/ape/yukiko/data/spctfilt/"
      @data =
        Ape.new("#{dir_path}#{$rezol}_#{$expID}_#{filt_hash[filter]}filt.nc")

      file_name = Array.new
      Dir.foreach($gr2ncfile_path) { |item|
        file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/
      }
      @data_filt = Ape.new(file_name[-1])

      # PR 
      if $rezol == "T159L48_eml"
        gphys = __gr_tr_tppn_open[true,-801..-401] 
        filt = @data.gphys_open("tr_tppn_#{filter}filt")[true,-801..-401] 
      else
        gphys = __gr_tr_tppn_open[true,-401..-1] 
        filt = @data.gphys_open("tr_tppn_#{filter}filt")[true,-401..-1] 
      end
      filt = filt.rename("tr_tppn_#{filter}filt_op")
      gphys = gphys.rename("tr_tppn_#{filter}filt_op")
      mkfig_plot_comp(gphys.lon_lotate, filt.lon_lotate)
      
      # T (0.57)
      if $rezol == "T159L48_eml"
        gphys = __gr_tr_lat0_open("T",0.57)[true,-801..-401] 
        filt = 
          @data.gphys_open("tr_t_#{filter}filt").
          cut(true,0.57,true)[true,-801..-401] 
      else
        gphys = __gr_tr_lat0_open("T",0.57)[true,-401..-1] 
        filt = 
          @data.gphys_open("tr_t_#{filter}filt").
          cut(true,0.57,true)[true,-401..-1] 
      end
      gphys = gphys.rename("tr_t570s_#{filter}filt_op")
      filt = filt.rename("tr_t570s_#{filter}filt_op")
      gphys = (gphys - gphys.mean(0)).
        add_lost_axes("(diff) from (mean) zonal")
      filt = (filt - filt.mean(0))
      mkfig_plot_comp(gphys.lon_lotate, filt.lon_lotate)

      # Q (0.85)
      if $rezol == "T159L48_eml"
        gphys = __gr_tr_lat0_open("Q",0.85)[true,-801..-401] 
        filt = 
          @data.gphys_open("tr_q_#{filter}filt").
          cut(true,0.85,true)[true,-801..-401] 
      else
        gphys = __gr_tr_lat0_open("Q",0.85)[true,-401..-1] 
        filt = 
          @data.gphys_open("tr_q_#{filter}filt").
          cut(true,0.85,true)[true,-401..-1] 
      end
      filt = filt.rename("tr_q850s_#{filter}filt_op")
      gphys = gphys.rename("tr_q850s_#{filter}filt_op")        
      gphys = (gphys - gphys.mean(0)).
        add_lost_axes("(diff) from (mean) zonal")
      filt = (filt - filt.mean(0))
      mkfig_plot_comp(gphys.lon_lotate, filt.lon_lotate)
      
      # U (0.25)
      if $rezol == "T159L48_eml"
        gphys = __gr_tr_lat0_open("U",0.25)[true,-801..-401] 
        filt = 
          @data.gphys_open("tr_u_#{filter}filt").
          cut(true,0.25,true)[true,-801..-401] 
      else
        gphys = __gr_tr_lat0_open("U",0.25)[true,-401..-1] 
        filt = 
          @data.gphys_open("tr_u_#{filter}filt").
          cut(true,0.25,true)[true,-401..-1] 
      end
      filt = filt.rename("tr_u250s_#{filter}filt_op")
      gphys = gphys.rename("tr_u250s_#{filter}filt_op")
      gphys = (gphys - gphys.mean(0)).
        add_lost_axes("(diff) from (mean) zonal")
      filt = (filt - filt.mean(0))
      mkfig_plot_comp(gphys.lon_lotate, filt.lon_lotate)
      
    }    

  end

  def __gr_tr_tppn_open
    g = @data_filt.gphys_open("tr_tppn").cut(true,0,true)
    g = g.set_lost_axes(g.lost_axes.to_s.sub("y=","lat=") )
    gphys = __gr2_gphys_make(g,0)
    return gphys
  end

  def __gr_tr_lat0_open(var,sigma)
    g = @data_filt.gphys_open(var).cut(true,sigma,true)
    lost_axis = [
      g.data.get_att("lost_axis").sub("y=","lat="),
      g.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")]
    g = g.set_lost_axes(lost_axis)
    
    gphys = __gr2_gphys_make(g,sigma)
    return gphys
  end
end