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

$local_path = '/home/yukiko/work/ape/yukiko/lib'
$: << $local_path
$fig_path = "/home/yukiko/tr_request/figs/tmp/"

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

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


def set_directry
  $groupid = $groupid_hash[$rezol]
  $ncfile_path = "/home/yukiko/tr_request/#{$rezol}/"
  $file_label = "#{$rezol}_#{$expID}"
end


END{

a = ape_mkfig_new(3)

#  ["AGUforAPE", "ECMWF", "GSFC", "LASG", "NCAR" ].each{ |resol|
#  ["ECMWF_07"].each{ |resol|
  ["CSIRO","ECMWF_07"].each{ |resol|

    sstid = ["control"]
    sstid.each{ |item|
      $expID = item
      $rezol = resol
      set_directry
      
##      a.nc_tr00
      a.nc_tr_qtuz

    }
  }

}

class Ape_mkfig

  # #{$groupid}_TR_control.nc
  def nc_tr00

    @data = Ape.new("#{$ncfile_path}#{$groupid}_TR00_#{$expID}.nc")
    @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR00_#{$expID}.nc")

#    @data = Ape.new("#{$ncfile_path}#{$groupid}_TR08_#{$expID}.nc")
#    @data_tmp = Ape.new("#{$ncfile_path}#{$groupid}_TR08_#{$expID}.nc")

    @data_tmp.netcdf_open.var_names.each { |name|
      
      if name =~ /lw_toa|mslp|500|250|ps|slh|ssh/
          if name =~ /lw_toa|mslp|ps/
            lost_axis = []
          elsif name =~ /250/
            lost_axis = ["plev=25000 Pa"]
          elsif name =~ /500/
            lost_axis = ["plev=50000 Pa"]
          end
          
          gphys = @data.gphys_open(name)
          
	  gphys =  gphys.cut(true,0,true)[true,-401..-1]
	  grid_1 =
	    Axis.new().
	    set_pos(VArray.new(NArray.sfloat(401).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 ).set_lost_axes(gphys.lost_axes)
  
         if gphys.name == "tr_ps" 
           gphys = gphys.rename("tr_mslp") 
           name = "tr_mslp"
         end

        if $rezol == "CSIRO"
           if name == "tr_lw_toa"
             gphys = gphys.set_att("ape_name","toa_outgoing_longwave_flux").
               set_att("units","W m-2")
           elsif name == "tr_mslp"
             gphys = gphys.set_att("ape_name","air_pressure_at_sea_level").
               set_att("units","Pa")
           elsif name == "tr_slh"
             gphys = gphys.set_att("ape_name","surface_upward_latent_heat_flux").
               set_att("units","W m-2")
           elsif name == "tr_ssh"
             gphys = gphys.set_att("ape_name","surface_upward_sensible_heat_flux").
               set_att("units","W m-2")
           end
        end

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

	  # 時空間スペクトル
	  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..80])

        elsif name =~ /tppn/

          gphys = @data.gphys_open(name)
  
	  # x-t ダイヤグラム
	    gphys = gphys.cut(true,0,true)[true,-401..-1]
     	    grid_1 =
	      Axis.new().
	      set_pos(VArray.new(NArray.sfloat(401).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 ).set_lost_axes(gphys.lost_axes)
	  
#	  if $rezol == "NCAR"
#	    gphys = gphys * 1000
#	  end
        if $rezol == "NCAR"
          gphys = gphys.set_att("ape_name", 
                                gphys.data.get_att("ape_name").
                                  gsub("kg m-2 s-1","precipitation flux")  )
        elseif $rezol == "CSIRO"
          gphys = gphys.set_att("ape_name","precipitation flux").
            set_att("units","kg m-2 s-1")
        end

  	  mkfig_plot(gphys.lon_lotate)
#	  mkfig_plot(gphys.lon_lotate.rename("tr_tppn_mono").set_lost_axes(""))


	  # 時空間スペクトル
	  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..80])
        end
    }

  end

  # #{$groupid}_TR_control.nc
  def nc_tr_qtuz
    
    name_array = [
      [13,"tr_t",25000,"tr_t250","air_temperature","K"], 
      [13,"tr_t",15000,"tr_t150","air_temperature","K"], 
      [13,"tr_t",60000,"tr_t600","air_temperature","K"], 
#---
      [16,"tr_q",92500,"tr_q925","specific_humidity","1"], 
      [16,"tr_q",85000,"tr_q850","specific_humidity","1"], 
      [16,"tr_q",70000,"tr_q700","specific_humidity","1"], 
#---
      [11,"tr_u",25000,"tr_u250","eastward_wind","m s-1"], 
      [11,"tr_u",15000,"tr_u150","eastward_wind","m s-1"], 
      [11,"tr_u",85000,"tr_u850","eastward_wind","m s-1"], 
      [11,"tr_u",60000,"tr_u600","eastward_wind","m s-1"], 
      [11,"tr_u",92500,"tr_u925","eastward_wind","m s-1"], 
#---
      [12,"tr_v",25000,"tr_v250","northward_wind","m s-1"], 
      [12,"tr_v",15000,"tr_v150","northward_wind","m s-1"], 
      [12,"tr_v",85000,"tr_v850","northward_wind","m s-1"], 
      [12,"tr_v",60000,"tr_v600","northward_wind","m s-1"], 
      [12,"tr_v",92500,"tr_v925","northward_wind","m s-1"], 
#---
      [15,"tr_z",25000,"tr_phi250","geopotential_height","m"], 
      [15,"tr_z",15000,"tr_phi150","geopotential_height","m"],
      [15,"tr_z",85000,"tr_phi850","geopotential_height","m"]
    ]
    
    name_array.length.times{ |num| 

      @data = Ape.new("#{$ncfile_path}#{$groupid}_TR#{name_array[num][0]}_#{$expID}.nc")


    lev = name_array[num][2]
    lev = lev /100    if $rezol =~ /ECMWF/ || $rezol == "LASG"|| $rezol == "CSIRO"
    opn_name = name_array[num][1]
    opn_name = "tr_phi" if $rezol == "AGUforAPE" && opn_name == "tr_z"

    gphys = @data.gphys_open(opn_name)
    gphys =  gphys.cut(true,0,lev,true)[true,-401..-1]

    gphys = gphys.rename(name_array[num][3]).
	    set_att("ape_name", name_array[num][4]).
	    set_att("units", name_array[num][5])

    grid_1 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(401).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 ).set_lost_axes(gphys.lost_axes)
  
  gphys = gphys /9.79764 if $rezol == "NCAR" && opn_name == "tr_z"

    mkfig_plot((gphys - gphys.mean(0)).
                 add_lost_axes("plev=#{name_array[num][2]} Pa").
                 add_lost_axes("(diff) from (mean) zonal").lon_lotate)
  
    # 時空間スペクトル
    gphys = gphys.stspct_fft("#{name_array[num][3]}_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("plev=#{name_array[num][2]} Pa")
	  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


end