#!/usr/bin/env ruby

=begin

表題: AGUforAPE NetCDF submitted PCMDI のお絵描き

履歴: 2004/02/13 yukiko@ep.sci.hokudai.ac.jp

irb_command: load "~/tmp/ape-data/lib/ape-view.rb"
irbrc: ~/.irbrc

=end

# ----------------------------------------------
# ps, gif ファイルの格納場所

# $rezol = "T39L48_eml"
$rezol = "T39L96_eml"

id = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
# $expID = "3keq"  # $expID = id[4]
$expID = id[0]


$ncfile_path = "/home/yukiko/tmp/ape-data/NetCDF/#{$rezol}/"
$file_label = "#{$rezol}_#{$expID}"
$fig_path = "/home/yukiko/tmp/ape-data/figs/tmp/"

$grfile_path = "/home/yukiko/tmp/ape-data/GrADS/#{$rezol}_expID01/"
$gr2ncfile_path = "/home/yukiko/tmp/ape-data/yukiko/#{$rezol}_expID01/"

# $ncfile_path = "/mnt/NetCDF/#{$rezol}/"
# $file_label = "T39L48_eml_control"

$grfile_path = "/work11/ape/GrADS/#{$rezol}_expID01/"
$gr2ncfile_path = "/work11/ape/yukiko/data/#{$rezol}_expID01/"

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

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

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

load "#{$local_path}/lib-ape-view.rb"
# require "lib-ape-view"
require "libgphys-e"
require "/usr/lib/ruby/1.6/irb/completion.rb" # tab で関数補完

# ----------------------------------------------
# irb 用ローカルメソッド (そもそもローカルだけど)
# ややこしくなるから irb でのみ使用すること 


# 読み込み時のメッセージ
def init_put

  print  <<EOF 

  ==========================================================
   Hello!! APE (Aqua-Planet experiment project) mapping irb
   Copyright: AGU for APE 2004 (YAMADA Yukiko)
  ---------------------------------------------------------

   * load file   : #{$local_path}/ape-view.rb
   * command help: ape_help 

  ==========================================================

EOF

end

# コマンドヘルプ
def ape_help(num=0)
  if num==1
    Ape.help
  elsif num==2
    Ape_mkfig.help
  else
    print  <<EOF 

  ==========================================================
   irb 用基本コマンド
  ---------------------------------------------------------
   * ape_help(num=0)  : help 出力 
      * num=0 : このヘルプメッセージ 
      * num=1 : class Ape のコマンド (lib-ape-view.rb) 
      * num=2 : class Ape_mkfig のコマンド (ape-view.rb) 
   * ape_source       : ape-view.rb の再読み込み
   * ape_set          : 読み込む ncfile 等の設定表示
   * ape_ls           : \"ls ./\" 出力. $nc に ncfile 名を代入.  
   * ape_new(file="#{$ncfile_path}AGUforAPE-03a_GT_#{$expID}.nc")
     : class Ape の初期化. file は指定しないと ape_set 値を適当に読む(未実装)  
   * ape_mkfig_new(num=1) : class Ape_mkfig の初期化. num は出力先 
      * num=1 : X
      * num=2 : ps
      * num=3 : gif (dump して生成)
  ==========================================================

EOF
  end

end

# 読み込みファイル等の確認
def ape_set
  print  <<EOF

  ==========================================================
   解像度     : $rezol       = #{$rezol}
   実験ID     : $expID       = #{$expID}
   ラベル     : $file_label  = #{$file_label}
   ncfile格納 : $ncfile_path = #{$ncfile_path}
   gif,ps格納 : $fig_path    = #{$fig_path}
  ==========================================================

EOF
end


# ライブラリのロード (irb 立ち上げ時に一度読み込まれてる)
def ape_source
  load "#{$local_path}/ape-view.rb"
end

# $nc[Array] に nc ファイル名代入
def ape_ls
  Ape.file_name_get
end

# ape 基本クラス初期化
def ape_new(file="#{$ncfile_path}AGUforAPE-03a_GT_#{$expID}.nc")
  Ape.new(file)
end

# ape 応用クラス初期化
def ape_mkfig_new(num=1)
  Ape_mkfig.new(num)
end

# Ape Grads -- NetCDF コンバータクラス初期化
def ape_mknetcdf_new
  Ape_mknetcdf.new
end


# 素直に ncdump -h
def ape_ncdump(file)
  Ape.ncdump(file)
end

# Ape メソッドエイリアス
class Ape
  def p(name) 
    gphys = gphys_open(name)
    plot(gphys)
  end
  alias go gphys_open
  alias no netcdf_open
  alias l var_list
end

# ape 
alias help ape_help
alias source ape_source
alias ls ape_ls

init_put

# -----------------------------------------------
# Ape NetCDF お絵描き応用クラス

class Ape_mkfig

  # class method ------------------------------------------

  def Ape_mkfig.help
    print  <<EOF

  ==========================================================
   class Ape_mkfig のメソッド (ape-view.rb) 
  ---------------------------------------------------------
   * Ape_mkfig.new(num)
   * mkfig_plot(gphys)
   * mkfig_plot_all
   * nc_varlist(file)
   * nc_gt, nc_mf, nc_ml, nc_sh, nc_sh_zonal, nc_tr .....
  ==========================================================

EOF

  end

  # method -----------------------------------------------

  def initialize(num)
    @fig_num = num
  end

  def help
    Ape_mkfig.help
  end

  def nc_varlist(file)
    @data = Ape.new(file)
    @data.netcdf_open.var_names.each { |item|
      unless item =~ /lon|lat|plev|time/ 
	gphys = @data.gphys_open(item)
	var_name = gphys.data.get_att("ape_name").gsub("_", " ")
	print "#{gphys.name}: #{var_name}\n"
      end
    }
  end

  def mkfig_plot(gphys)
    if @fig_num == 3
      @data.mkdump(gphys)
    elsif @fig_num == 2
      @data.mkps(gphys)
    else
      @data.plot(gphys)
    end
  end

  def nc_plot_all
    $file_name = nil ;  $pre_file = nil

#    nc_sh
#    nc_sh
#    nc_tr
#    nc_gt
    nc_ml
    nc_ml
#    nc_mf
#    nc_sh_zonal
#    nc_sh_zonal

=begin

    nc_sh_anm
    nc_sh_anm
    nc_gt_anm
    nc_ml_anm
    nc_mf_anm
    nc_sh_zonal_anm

    nc_tr
    nc_sh
    nc_gt
    nc_ml
    nc_mf
    nc_sh_zonal
    nc_sh_zonal
=end

  end

  # AGUforAPE-03a_GT_control.nc
  def nc_gt
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_GT_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item| 
      unless item =~ "time"  || item =~ "gt_slwu" # gt_slwu は一定値
	gphys = @data.gphys_open(item)
	mkfig_plot(gphys)
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("gt_cppn")
    g2 = @data.gphys_open("gt_dppn")
    gphys = (g1 + g2).rename("gt_tppn").
      set_att("ape_name", "total_precipitation_flux")
    mkfig_plot(gphys)
  end

  # AGUforAPE-03a_ML_control.nc
  def nc_ml
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_ML_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "lon" || item =~ "plev"  
	gphys = @data.gphys_open(item)
	mkfig_plot(gphys.mean(0))
      end
    }
    # mass_stream_function
    gphys = @data.gphys_open("ml_v")
    gphys = gphys.mean(0).strm_function("ml_psi").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1") 
    mkfig_plot(gphys)
  end

  # AGUforAPE-03a_TR_control.nc
  def nc_tr
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_TR_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |name|
      if name =~ /lw_toa|mslp|500|250/
	if name =~ /lw_toa|mslp/
	  lost_axis = []
	elsif name =~ /250/
	  lost_axis = ["plev=25000 Pa"]
	elsif name =~ /500/
	  lost_axis = ["plev=50000 Pa"]
	end

	# x-t ダイヤグラム
	gphys = @data.gphys_open(name)
	gphys =  gphys.cut(true,0,true)[true,-400..-1]
	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")})2 day").
	  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)/4-1)..((dim1+1)/4*3-1),0..(dim2/4)])

      elsif name =~ /tppn/

	# x-t ダイヤグラム
	gphys = @data.gphys_open("tr_tppn")
	gphys = gphys.cut(true,0,true)[true,-400..-1]
	mkfig_plot(gphys.lon_lotate)

	# 時空間スペクトル
	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")})2 day").
	  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)/4-1)..((dim1+1)/4*3-1),0..(dim2/4)])
      end
    }

  end

  #  AGUforAPE-03a_SH_control.nc
  def nc_sh
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_SH_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys = @data.gphys_open(item)
	gphys = gphys.lon_lotate
	mkfig_plot(gphys)
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2).rename("sh_tppn").
      set_att("ape_name", "total_precipitation_flux")
    gphys = gphys.lon_lotate
    mkfig_plot(gphys)
  end

  #  AGUforAPE-03a_SH_control.nc@zonal_mean
  def nc_sh_zonal
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_SH_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys = @data.gphys_open(item)
	gphys = gphys.rename("#{gphys.name}_zonal")
	mkfig_plot(gphys.mean(0))
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2).rename("sh_tppn_zonal").
      set_att("ape_name", "total_precipitation_flux")
    mkfig_plot(gphys.mean(0))
  end
  
  #  AGUforAPE-03a_MF_control.nc
  def nc_mf
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_MF_#{$expID}.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "plev" 
	gphys = @data.gphys_open(item)
	mkfig_plot(gphys)
      end
    }
  end

  # anomary -----------------------------------------------------

  def nc_gt_anm
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_GT_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}AGUforAPE-03a_GT_control.nc")
    @data.netcdf_open.var_names.each { |item| 
      unless item =~ "time"  || item =~ "gt_slwu" # gt_slwu は一定値
	gphys_con = @data_con.gphys_open(item).mean(0)
	lost_axis = ["(diff) from (mean) time of control"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys - gphys_con)
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("gt_cppn")
    g2 = @data.gphys_open("gt_dppn")
    gphys = (g1 + g2)
    g1 = @data_con.gphys_open("gt_cppn")
    g2 = @data_con.gphys_open("gt_dppn")
    gphys_con = (g1 + g2).mean(0)
    lost_axis = ["(diff) from (mean) time of control"]
    gphys = (gphys - gphys_con).rename("gt_tppn_anm").
      set_att("ape_name", "total_precipitation_flux").
      add_lost_axes(lost_axis)
    mkfig_plot(gphys)
  end

  def nc_ml_anm
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_ML_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}AGUforAPE-03a_ML_control.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "lon" || item =~ "plev"  
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) control"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys.mean(0) - gphys_con.mean(0))
      end
    }
    # mass_stream_function
    lost_axis = ["(diff) from (mean) control"]
    gphys = @data.gphys_open("ml_v")
    gphys = gphys.add_lost_axes(lost_axis).mean(0).strm_function("ml_psi_anm").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1")
    gphys_con = @data_con.gphys_open("ml_v")
    gphys_con = gphys_con.mean(0).strm_function("ml_psi_anm").
      set_att("ape_name","mass_stream_function").set_att("units","Kg s-1") 
    mkfig_plot(gphys-gphys_con)
  end

  def nc_sh_anm
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_SH_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}AGUforAPE-03a_SH_control.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) zonal of control"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = (gphys - gphys_con.mean(0))
	gphys = gphys.lon_lotate.rename("#{gphys.name}_anm")
	mkfig_plot(gphys)
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2)
    g1 = @data_con.gphys_open("sh_cppn")
    g2 = @data_con.gphys_open("sh_dppn")
    gphys_con = (g1 + g2)
    gphys = (gphys - gphys_con.mean(0))
    lost_axis = ["(diff) from (mean) zonal of control"]
    gphys = gphys.lon_lotate.rename("sh_tppn_anm").
      set_att("ape_name", "total_precipitation_flux").
      add_lost_axes(lost_axis)
    mkfig_plot(gphys)
  end

  def nc_sh_zonal_anm
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_SH_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}AGUforAPE-03a_SH_control.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lon" || item =~ "lat" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) control"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_zonal_anm")
	mkfig_plot(gphys.mean(0)- gphys_con.mean(0))
      end
    }
    # total_precipitation_flux
    g1 = @data.gphys_open("sh_cppn")
    g2 = @data.gphys_open("sh_dppn")
    gphys = (g1 + g2)
    g1 = @data_con.gphys_open("sh_cppn")
    g2 = @data_con.gphys_open("sh_dppn")
    gphys_con = (g1 + g2)
    gphys = (gphys.mean(0) - gphys_con.mean(0))
    lost_axis = ["(diff) from (mean) control"]
    gphys = gphys.rename("sh_tppn_zonal_anm").
      set_att("ape_name", "total_precipitation_flux").
      add_lost_axes(lost_axis)
    mkfig_plot(gphys)
  end
  
  def nc_mf_anm
    @data = Ape.new("#{$ncfile_path}AGUforAPE-03a_MF_#{$expID}.nc")
    @data_con = Ape.new("#{$ncfile_path}AGUforAPE-03a_MF_control.nc")
    @data.netcdf_open.var_names.each { |item|
      unless item =~ "lat" || item =~ "plev" 
	gphys_con = @data_con.gphys_open(item)
	lost_axis = ["(diff) from (mean) control"]
	gphys = @data.gphys_open(item).add_lost_axes(lost_axis)
	gphys = gphys.rename("#{gphys.name}_anm")
	mkfig_plot(gphys - gphys_con)
      end
    }
  end

  # GrADS -----------------------------------------------------
  
  def gr_tr
    @data = Ape.new("PRCPC.ctl")
    gphys = @data.gphys_open("PRCPC").cut(true,0,1,true)[true,-400..-1]
    gphys_add = Ape.new("PRCPL.ctl").gphys_open("PRCPL").
      cut(true,0,1,true)[true,-400..-1]
    gphys = gphys + gphys_add
    lost_axes = Array.new
    gphys.lost_axes.each{ |i| 
      lost_axes.push( i.sub("y=","lat=") ) unless i =~ /z=/ }

    # x-t ダイヤグラム
    gphys = gphys.rename("tr_tppn").
      set_att("ape_name","precipitation flux").
      set_att("units","kg m-2 s-1").
      set_lost_axes(lost_axes)
    mkfig_plot(gphys.lon_lotate)

    # 時空間スペクトル
    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")})2 day").
      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)/4-1)..((dim1+1)/4*3-1),0..(dim2/4)])

  end

  def gr_test
    @data = Ape.new("PRCPC.ctl")
    gphys = @data.gphys_open("PRCPC").cut(true,0,1,true)[true,-400..-1]
    gphys_add = Ape.new("PRCPL.ctl").gphys_open("PRCPL").
      cut(true,0,1,true)[true,-400..-1]
    gphys = gphys + gphys_add
    lost_axes = Array.new
    gphys.lost_axes.each{ |i| 
      lost_axes.push( i.sub("y=","lat=") ) unless i =~ /z=/ }

    gphys = gphys.rename("tr_tppn").
      set_att("ape_name","precipitation flux").
      set_att("units","kg m-2 s-1")


    # nc に一度落す
    file = NetCDF.create('tmp.nc')
    GPhys::NetCDF_IO.write(file, gphys)
    file.close

    # x-t ダイヤグラム
    @data = Ape.new("tmp.nc")
    gphys = @data.gphys_open("tr_tppn").set_lost_axes(lost_axes)
    mkfig_plot(gphys.lon_lotate)

    # 時空間スペクトル
    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")})2 day").
      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)/4-1)..((dim1+1)/4*3-1),0..(dim2/4)])

  end

  def gr_test1

    ["0013","0014"].each{ |num|

      @data = Ape.new("#{num}/PRCPC.ctl")
      gphys = @data.gphys_open("PRCPC").cut(true,0,1,true)
      gphys_add = Ape.new("#{num}/PRCPL.ctl").gphys_open("PRCPL").
	cut(true,0,1,true)
      gphys = gphys + gphys_add
      lost_axes = Array.new
      gphys.lost_axes.each{ |i| 
	lost_axes.push( i.sub("y=","lat=") ) unless i =~ /z=/ }
      
      gphys = gphys.rename("tr_tppn").
	set_att("ape_name","precipitation flux").
	set_att("units","kg m-2 s-1"). 
	set_att("lost_axis",lost_axes.to_s)

      # nc に一度落す
      file = NetCDF.create("tmp_#{num}.nc")
      GPhys::NetCDF_IO.write(file, gphys)
      file.close

    }
    
    # x-t ダイヤグラム
    @data = Ape.new(["tmp_0013.nc","tmp_0014.nc"])
    gphys = GPhys::NetCDF_IO.open(["tmp_0013.nc","tmp_0014.nc"], "tr_tppn")
    gphys = gphys[true,-400..-1].set_lost_axes(gphys.data.get_att("lost_axis"))
    mkfig_plot(gphys.lon_lotate)

    # 時空間スペクトル
    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")})2 day").
      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)/4-1)..((dim1+1)/4*3-1),0..(dim2/4)])

  end

end


# -----------------------------------------------
# Ape Grads -- NetCDF コンバータ

class Ape_mknetcdf

  def mknetcdf
  
    grdir = Array.new
    Dir.foreach($grfile_path) { |item|
      if item =~ /^00/ 
	grdir.push(item)
	print "#{item}\n"
      end
    }
    
    # [PRCPC, PRCPL, PS]
    grdir.each{ |num|
      
      file = NetCDF.create("#{$gr2ncfile_path}/#{$rezol}_expID01_#{num}.nc")
      
      @data = Ape.new("#{$grfile_path}/#{num}/PRCPC.ctl")
      gphys = @data.gphys_open("PRCPC").cut(true,0,1,true)
      gphys_add = Ape.new("#{$grfile_path}/#{num}/PRCPL.ctl").
	gphys_open("PRCPL").cut(true,0,1,true)
      gphys = gphys + gphys_add
      lost_axes = Array.new
      gphys.lost_axes.each{ |i| 
	lost_axes.push( i.sub("y=","lat=") ) unless i =~ /z=/ }
      gphys = gphys.rename("tr_tppn").
	set_att("ape_name","precipitation flux").
	set_att("units","kg m-2 s-1"). 
	set_att("lost_axis",lost_axes.to_s)
      GPhys::NetCDF_IO.write(file, gphys)
      
      file.close
    }
    
  end

end