#!/usr/bin/env ruby

=begin

表題: AGUforAPE GrADS to NetCDF

履歴: 2005/02/18 yukiko＠ep.sci.hokudai.ac.jp

=end

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


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


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

END{

  set_opt

  print "$expid = #{$expid}, $resol = #{$resol}\n"

  ncfile_pf_make
  print "PF NetCDF make end\n"


}

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

def set_opt

  $resol = "T39L48_eml"
  $expid = 1

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

  expID = "expID0#{$expid}"
  id = ["control","peaked","flat","Qobs","control-5N","1keq","3keq","3kw1"]
  $sst = id[$expid.to_i - 1]

#  $ncdir = "/home/yukiko/work/ape/NetCDF/"
  $ncdir = "/home/yukiko/work/ape/NetCDF/#{$resol}/"
#  $ncdir = "/work11/ape/NetCDF/#{$resol}/"
  $ncfile_head = "AGUforAPE-03a"
  $grddir = "/home/yukiko/work/ape/GrADS/#{$resol}_#{expID}/mean/"
#  $grddir = "/mnt/hgfs/scsi/#{$resol}_#{expID}/mean/"
  puts $grddir

end


# sub() ========================================================
def ncfile_pf_make

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

  item = "DTRADS"
  gphys = GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl", item)

  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 = g1.set_att("standard_name","latitude")
  g1 = g1.set_att("units","degrees_north")
  g1 = Axis.new().set_pos(g1)

  g2 = VArray.new(gphys.grid_copy.axis(2).pos.val ).rename("lev")
  g2 = g2.set_att("units","sigma_level")
  g2 = g2.set_att("standard_name","atmosphere_simga_coordinate")
  g2 = Axis.new().set_pos(g2)

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


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

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

  $cell_methods = "time: mean (accumulation at every time-step comment: averaged over 3-year period following 6-month spin-up)"

  gphys_ary = []

  # sh_ps ----------------------------------------------------
  
  print "  sh_ps start #{Time.now}\n"
  item = 'PS'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)[true,true,0]

  data = VArray.new(gphys.data.val).rename("sh_ps")
  data = data.set_att("units","Pa")
  data = data.set_att("ape_name","surface_air_pressure")
  data = data.set_att("cell_methods", "time: mean (interval: 6 hours comment: averaged over 3-year period following 6-month spin-up)")
  gphys_ary.push(GPhys.new(grid_ps,data))

  # pf_t ----------------------------------------------------
  
  print "  pf_t start #{Time.now}\n"
#  item_ary = ['DTRADS','DTRADL','DTVDF','DTCUM','DTLSC','DTDYN','DTDCA']
  item_ary = ['DTRADS','DTRADL','DTVDF','DTCUM','DTLSC','DTDCA']

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item_ary[0]}.ctl",item_ary[0]).cut(true,true,true,0) * 0.0

  item_ary.each{ |item|
    work = 
    	GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
	cut(true,true,true,0)
    gphys = gphys + work.copy
  }

  data = VArray.new(gphys.data.val).rename("pf_t")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Total")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_t_sw ----------------------------------------------------
  
  print "  pf_t_sw start #{Time.now}\n"
  item = 'DTRADS'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_t_sw")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Short wave")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_t_lw ----------------------------------------------------
  
  print "  pf_t_lw start #{Time.now}\n"
  item = 'DTRADL'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_t_lw")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Long wave")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_t_turb ----------------------------------------------------
  
  print "  pf_t_turb start #{Time.now}\n"
  item = 'DTVDF'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_t_turb")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Turbulence")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_t_conv ----------------------------------------------------
  
  print "  pf_t_conv start #{Time.now}\n"
  item = 'DTCUM'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_t_conv")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Convection")
  data = data.set_att("long_name","Temperature tendency - Cumulus Convection")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_t_cld ----------------------------------------------------
  
  print "  pf_t_cld start #{Time.now}\n"
  item = 'DTLSC'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_t_cld")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Cloud")
  data = data.set_att("long_name","Temperature tendency - Large Scale Condensation")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

=begin
  # pf_t_dca ----------------------------------------------------
  
  print "  pf_t_dca start #{Time.now}\n"
  item = 'DTDCA'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_t_dca")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Dry Convection Adjustment")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_t_dyn ----------------------------------------------------
  
  print "  pf_t_dyn start #{Time.now}\n"
  item = 'DTDYN'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_t_dyn")
  data = data.set_att("units","K s-1")
  data = data.set_att("ape_name","Temperature tendency - Dynamics")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))
=end

  # pf_q ----------------------------------------------------
  
  print "  pf_q start #{Time.now}\n"
#  item_ary = ['DQVDF', 'DQCUM', 'DQLSC', 'SRMNQ', 'DQDYN', 'DQDCA']
  item_ary = ['DQVDF', 'DQCUM', 'DQLSC', 'SRMNQ', 'DQDCA']

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item_ary[0]}.ctl",item_ary[0]).cut(true,true,true,0) * 0.0

  item_ary.each{ |item|
    work = 
    	GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
	cut(true,true,true,0)
    gphys = gphys + work.copy
  }

  data = VArray.new(gphys.data.val).rename("pf_q")
  data = data.set_att("units","kg kg-1 s-1")
  data = data.set_att("ape_name","Absolute humidity tendency - Total")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_q_turb ----------------------------------------------------
  
  print "  pf_q_turb start #{Time.now}\n"
  item = 'DQVDF'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_q_turb")
  data = data.set_att("units","kg kg-1 s-1")
  data = data.set_att("ape_name","Absolute humidity tendency - Turbulence")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_q_conv ----------------------------------------------------
  
  print "  pf_q_conv start #{Time.now}\n"
  item = 'DQCUM'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_q_conv")
  data = data.set_att("units","kg kg-1 s-1")
  data = data.set_att("ape_name","Absolute humidity tendency - Convection")
  data = data.set_att("long_name","Absolute humidity tendency - Cumulus Convection")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_q_cld ----------------------------------------------------
  
  print "  pf_q_cld start #{Time.now}\n"
  item = 'DQLSC'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_q_cld")
  data = data.set_att("units","kg kg-1 s-1")
  data = data.set_att("ape_name","Absolute humidity tendency - Cloud")
  data = data.set_att("long_name","Absolute humidity tendency - Large Scale Condensation")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_q_negq ----------------------------------------------------
  
  print "  pf_q_negq start #{Time.now}\n"
  item = 'SRMNQ'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_q_negq")
  data = data.set_att("units","kg kg-1 s-1")
  data = data.set_att("ape_name","Absolute humidity tendency - Negative q fix")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

=begin
  # pf_q_dca ----------------------------------------------------
  
  print "  pf_q_dca start #{Time.now}\n"
  item = 'DQDCA'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_q_dca")
  data = data.set_att("units","kg kg-1 s-1")
  data = data.set_att("long_name","Absolute humidity tendency - Dry Convection Adjustment")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_q_dyn ----------------------------------------------------
  
  print "  pf_q_dyn start #{Time.now}\n"
  item = 'DQDYN'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_q_dyn")
  data = data.set_att("units","kg kg-1 s-1")
  data = data.set_att("long_name","Absolute humidity tendency - Dynamics")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))
=end

  # pf_u ----------------------------------------------------
  
  print "  pf_u start #{Time.now}\n"
#  item_ary = ['DUVDF', 'DUGWD', 'DUDYN']
#  item_ary = ['DUVDF', 'DUGWD', 'DUDYN','DUCUM']
#  item_ary = ['DUPHY', 'DUDYN']
  item_ary = ['DUVDF', 'DUGWD', 'DUCUM']

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item_ary[0]}.ctl",item_ary[0]).cut(true,true,true,0) * 0.0

  item_ary.each{ |item|
    work = 
    	GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
	cut(true,true,true,0)
    gphys = gphys + work.copy
  }

  data = VArray.new(gphys.data.val).rename("pf_u")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Zonal velocity tendency - Total")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_u_turb ----------------------------------------------------
  
  print "  pf_u_turb start #{Time.now}\n"
  item = 'DUVDF'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_u_turb")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Zonal velocity tendency - Turbulence")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_u_conv ----------------------------------------------------
  
  print "  pf_u_conv start #{Time.now}\n"
  item = 'DUCUM'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_u_conv")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Zonal velocity tendency - Convection")
  data = data.set_att("long_name","Zonal velocity tendency - Cumulus Convection")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_u_gwd ----------------------------------------------------
  
  print "  pf_u_gwd start #{Time.now}\n"
  item = 'DUGWD'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_u_gwd")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Zonal velocity tendency - Gravity-wave drag")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

=begin
  # pf_u_dyn ----------------------------------------------------
  
  print "  pf_u_dyn start #{Time.now}\n"
  item = 'DUDYN'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_u_dyn")
  data = data.set_att("units","m s-2")
  data = data.set_att("long_name","Zonal velocity tendency -  Dynamics")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))
=end

  # pf_v ----------------------------------------------------
  
  print "  pf_v start #{Time.now}\n"
#  item_ary = ['DVVDF', 'DVGWD', 'DVDYN']
#  item_ary = ['DVVDF', 'DVGWD', 'DVDYN','DVCUM']
  item_ary = ['DVVDF', 'DVGWD', 'DVCUM']

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item_ary[0]}.ctl",item_ary[0]).cut(true,true,true,0) * 0.0

  item_ary.each{ |item|
    work = 
    	GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
	cut(true,true,true,0)
    gphys = gphys + work.copy
  }

  data = VArray.new(gphys.data.val).rename("pf_v")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Meridional velocity tendency - Total")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_v_turb ----------------------------------------------------
  
  print "  pf_v_turb start #{Time.now}\n"
  item = 'DVVDF'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_v_turb")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Meridional velocity tendency - Turbulence")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_v_conv ----------------------------------------------------
  
  print "  pf_v_conv start #{Time.now}\n"
  item = 'DVCUM'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_v_conv")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Meridional velocity tendency - Convection")
  data = data.set_att("long_name","Meridional velocity tendency - Cumulus Convection")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

  # pf_v_gwd ----------------------------------------------------
  
  print "  pf_v_gwd start #{Time.now}\n"
  item = 'DVGWD'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_v_gwd")
  data = data.set_att("units","m s-2")
  data = data.set_att("ape_name","Meridional velocity tendency - Gravity-wave drag")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))

=begin
  # pf_v_dyn ----------------------------------------------------
  
  print "  pf_v_dyn start #{Time.now}\n"
  item = 'DVDYN'

  gphys = 
    GPhys::GrADS_IO.open("#{$grddir}/mean#{item}.ctl",item).
    cut(true,true,true,0)

  data = VArray.new(gphys.data.val).rename("pf_v_dyn")
  data = data.set_att("units","m s-2")
  data = data.set_att("long_name","Meridional velocity tendency -  Dynamics")
  data = data.set_att("cell_methods", $cell_methods)
  gphys_ary.push(GPhys.new(grid,data))
=end

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

  print "  netcdf write start #{Time.now}\n"
  f=NetCDF.create("#{$ncdir}/#{$ncfile_head}_PF_#{$sst}.nc")

  londim = f.def_dim("lon",g0.pos.val.size)
  lon    = f.def_var("lon","float",[londim])
  lon.put_att("standard_name","longitude")
  lon.put_att("units","degrees_east")

  latdim = f.def_dim("lat",g1.pos.val.size)
  lat    = f.def_var("lat","float",[latdim])
  lat.put_att("standard_name","latitude")
  lat.put_att("units","degrees_north")

  levdim = f.def_dim("lev",g2.pos.val.size)
  lev    = f.def_var("lev","float",[levdim])
  lev.put_att("standard_name","atmosphere_simga_coordinate")
  lev.put_att("units","1")

  ilevdim = f.def_dim("ilev",$half_sigma.size)
  ilev    = f.def_var("ilev","float",[ilevdim])
  ilev.put_att("standard_name","atmosphere_simga_coordinate_half_level")
  ilev.put_att("units","1")

  f.enddef
 
  lon.put(g0.pos.val)
  lat.put(g1.pos.val)
  lev.put(g2.pos.val)
  ilev.put(NArray.to_na($half_sigma.reverse))

  f.redef

  gphys_ary.each{ |gphys|

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


  f.put_att("Conventions", "CF-1.0")
  f.put_att("title","Aqua Planet: Parameterization Forcing from  \'#{$sst}\' Experiment")
  f.put_att("history", "Original data produced: 2005/08/03 on ES; time-mean calculated for APE PF standard output")
  f.put_att("institution", "AGU for APE; AFES Working Team, GFD Dennou Club, University of Tokyo")
  f.put_att("source","AFES v1.15")
  f.close


end


$half_sigma = [ 
    0.0000,      0.29858E-02, 0.61058E-02, 0.93660E-02, 0.12773E-01,
    0.16332E-01, 0.20052E-01, 0.23939E-01, 0.28000E-01, 0.32281E-01,
    0.37217E-01, 0.42907E-01, 0.49468E-01, 0.57031E-01, 0.65751E-01,
    0.75804E-01, 0.87394E-01, 0.10076,     0.11616,     0.13392,
    0.15440,     0.17800,     0.20522,     0.23660,     0.27277,
    0.31448,     0.36256,     0.41799,     0.47952,     0.53869,
    0.59452,     0.64634,     0.69377,     0.73669,     0.77513,
    0.80927,     0.83938,     0.86577,     0.88878,     0.90877,
    0.92605,     0.94095,     0.95377,     0.96478,     0.97420,
    0.98225,     0.98913,     0.99500,     1.0000
 ]