#!/usr/bin/env ruby

=begin

表題: AGCM5 mradl TR NetCDF 作成ツール

履歴: 2005/04/18 yukiko@ep.sci.hokudai.ac.jp

=end

# 設定 =========================================================
# local load path
$local_path = '/work11/ape/yukiko/lib'
$: << $local_path


# 必要なライブラリ, モジュールの読み込み
require "lib-ape-base.rb"


# main() =======================================================

END{

#  $cum = "adj"
  $cum = "kuo"
  $expid = "mradlAd"
#  $expid = "con"

  $ncdir = "/work11/aqua3d/yukiko/mradl/data/test"
  $ncfile_head = "AGCM5-#{$cum}-#{$expid}"
  $gt3dir = "/work11/aqua3d/yukiko/mradl/data/test/#{$cum}700#{$expid}"

  ncfile_tr_make
  print "TR NetCDF make end\n"

}

# 設定 ========================================================

def set_opt


  $expid = 1

  # 実行オプション解析 
  $cum = "kuo"                             if ARGV.index("-kuo") 
  $cum = "kuo-nosc"                        if ARGV.index("-kuo-nosc") 
  $expid = ARGV[ARGV.index("-exp")+1]      if ARGV.index("-exp") 

  if $expid == "flat3keq" ||$expid == "Qobs3keq" || $expid =~ /H1998/ || $expid =~ /HS1986/
    expID = $expid
    $sst  = $expid
  else 
    expID = "expID0#{$expid}"
    id = ["control","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1"]
    $sst = id[$expid.to_i - 1]

  end

  $ncdir = "/work11/ape/NetCDF/AGCM5_#{$cum}/"
  $ncfile_head = "AGCM5-#{$cum}"
  $gt3dir = "/work11/ape/GTOOL3/#{$cum}_#{expID}"
#  $gt3dir = "/work11/ape/GTOOL3/AGCM5-#{$cum}_#{expID}"
#  $gt3dir = "/work11/ape/GTOOL3/AGCM5-#{$cum}_add_#{expID}"
  puts $gt3dir

end


# sub() ========================================================
def ncfile_tr_make

# grid =========================================================

  item = "ps_out.nc"
  gphys = GPhys::NetCDF_IO.open("#{$gt3dir}400/#{item}", "ps")

  g0 = VArray.new(gphys.grid_copy.axis(0).pos.val ).rename("lon")
  g0 = g0.set_att("standard_name","longitude")
  g0 = g0.set_att("units","degrees_east")
  g0 = Axis.new().set_pos(g0)

  g1 = VArray.new(gphys.grid_copy.axis(1).pos.val ).rename("lat")
  g1[-1..0] = g1[0..-1] # 緯度軸の反転
  g1 = g1.set_att("standard_name","latitude")
  g1 = g1.set_att("units","degrees_north")
  g1 = Axis.new().set_pos(g1)

  g2 = VArray.new( NArray.sfloat(1440).indgen!(1)*0.25 ).rename("time")
  g2 = g2.set_att("standard_name","time")
  g2 = g2.set_att("units","days since 0000-01-01")
  g2 = g2.set_att("bounds","time_bnds")
  g2 = Axis.new().set_pos(g2)

  grid = Grid.new(g0,g1,g2)

# variable ===================================================

  # variable common ---------------------------------------------

  $cell_methods = "time: mean (interval: accumulation at every time-step)"

  gphys_ary = []

  # rain ----------------------------------------------------
  
  print "  tr_tppn start\n"
  item = "rain_sum.nc"

  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "rain")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val*0.4*10**(-6)).rename("tr_tppn")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","kg m-2 s-1")
  data = data.set_att("ape_name","precipitation_flux")
  data = data.set_att("cell_methods", "time: mean (interval: accumulation at every time-step)")
  
  gphys_ary.push(GPhys.new(grid,data))
  
  # u_z13 ----------------------------------------------------
  
  print "  tr_u_z13 start\n"
  item = "u_out_z13.nc"

  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "u")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val).rename("tr_u_z13")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","m s-1")
  data = data.set_att("ape_name","eastward_wind")
  data = data.set_att("cell_methods", "time: point")
  
  gphys_ary.push(GPhys.new(grid,data))

  # u_z12 ----------------------------------------------------
  
  print "  tr_u_z12 start\n"
  item = "u_out_z12.nc"

  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "u")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val).rename("tr_u_z12")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","m s-1")
  data = data.set_att("ape_name","eastward_wind")
  data = data.set_att("cell_methods", "time: point")
  
  gphys_ary.push(GPhys.new(grid,data))

  
  # v_z13 ----------------------------------------------------
  
  print "  tr_v_z13 start\n"
  item = "v_out_z13.nc"

  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "v")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val).rename("tr_v_z13")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","m s-1")
  data = data.set_att("ape_name","northward_wind")
  data = data.set_att("cell_methods", "time: point")
  
  gphys_ary.push(GPhys.new(grid,data))

  # v_z12 ----------------------------------------------------
  
  print "  tr_v_z12 start\n"
  item = "v_out_z12.nc"

  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "v")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val).rename("tr_v_z12")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","m s-1")
  data = data.set_att("ape_name","northward_wind")
  data = data.set_att("cell_methods", "time: point")
  
  gphys_ary.push(GPhys.new(grid,data))
  
  # gph_z13 ----------------------------------------------------
  
  print "  tr_gph_z13 start\n"
  item = "gph_out_z13.nc"

  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "z")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val).rename("tr_gph_z13")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","m")
  data = data.set_att("ape_name","geopotential_height")
  data = data.set_att("cell_methods", "time: point")
  
  gphys_ary.push(GPhys.new(grid,data))
  
  # gph_z12 ----------------------------------------------------
  
  print "  tr_gph_z12 start\n"
  item = "gph_out_z12.nc"

  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "z")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val).rename("tr_gph_z12")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","m")
  data = data.set_att("ape_name","geopotential_height")
  data = data.set_att("cell_methods", "time: point")
  
  gphys_ary.push(GPhys.new(grid,data))
  
  # mslp ----------------------------------------------------
  
  print "  tr_mslp start\n"
  item = "ps_out.nc"
  
  file = ["#{$gt3dir}400/#{item}","#{$gt3dir}500/#{item}",
    "#{$gt3dir}600/#{item}","#{$gt3dir}700/#{item}"]
  gphys = GPhys::NetCDF_IO.open(file, "ps")[true,true,-1440..-1]
  
  data = VArray.new(gphys.data.val*100).rename("tr_mslp")
  data[true,-1..0,true] = data[true,0..-1,true] # 緯度軸の反転
  data = data.set_att("units","Pa")
  data = data.set_att("ape_name","air_pressure_at_sea_level")
  data = data.set_att("cell_methods", "time: point")
  
  gphys_ary.push(GPhys.new(grid,data))
  

  # NetCDF 書き込み ================================================

  print "  nc create start\n"
  f=NetCDF.create("#{$ncdir}/#{$ncfile_head}_TR.nc")

  gphys_ary.each{ |gphys|

    GPhys::NetCDF_IO.write( f, gphys )
    
  }

  f.put_att("Conventions", "CF-1.0")
  f.put_att("title","Aqua Planet: Transients from \'#{$expid}\' Experiment")
  f.put_att("history", "Original data produced: #{Time.now} on NIES SX6")
  f.put_att("institution", "GFD Dennou Club")
  f.put_att("source","agcm5.3")
  f.close

end