#!/usr/bin/env ruby

=begin

表題: AGUforAPE GrADS to NetCDF

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


=end

END{

# $rezol = "T39L48_eml"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T39L96_eml"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T39L24_eml"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T79L48_eml"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T39L48_non"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T159L48_eml"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T39L96_non"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T39L24_non"
# Ape_mknetcdf.new.mknetcdf

# $rezol = "T79L48_non"
# Ape_mknetcdf.new.mknetcdf

#$rezol = "T319L48_eml"
#Ape_mknetcdf.new.mknetcdf

# $rezol = "T159L48_non"
# $rezol = "T79L48_non"
# $rezol = "T39L24_non"
# Ape_mknetcdf.new.mknetcdf

 $rezol = "T39L48_ias"
 Ape_mknetcdf.new.mknetcdf

 $rezol = "T39L48_ksc"
 Ape_mknetcdf.new.mknetcdf

 $rezol = "T39L48_kuo"
 Ape_mknetcdf.new.mknetcdf

 $rezol = "T39L48_mca"
 Ape_mknetcdf.new.mknetcdf

 $rezol = "T39L48_non"
 Ape_mknetcdf.new.mknetcdf


}


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

id = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]

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

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

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

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


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

class Ape_mknetcdf

  def mknetcdf
  
    grfile_path = "/home/yukiko/work/ape/GrADS/#{$rezol}_expID01/"
    gr2ncfile_path = "/home/yukiko/work/ape/yukiko/data/#{$rezol}_expID01/"

    grdir = Array.new
    Dir.foreach(grfile_path) { |item|
      if item =~ /^00/ 
	grdir.push(item) unless item == "0001"
	print "#{item}\n"
      end
    }

    pwd_path = Dir.pwd

#    grdir = ["0014"]
    grdir.each{ |num|

      file = NetCDF.create("#{gr2ncfile_path}/#{$rezol}_expID01_horizontal_#{num}.nc")

      Dir.chdir("#{grfile_path}/#{num}")
      print Dir.pwd + "\n"

      # tr_tppn
      if  File.exist?("PRCPC.ctl") && File.exist?("PRCPC.ctl") then
	print "GrADS to NetCDF PRCPC & PRCPL\n"
	@data = Ape.new("#{grfile_path}/#{num}/PRCPC.ctl")
	gphys = @data.gphys_open("PRCPC").cut(true,-30..30,1,true)
	gphys_add = Ape.new("#{grfile_path}/#{num}/PRCPL.ctl").
	  gphys_open("PRCPL").cut(true,-30..30,1,true)
	gphys = gphys + gphys_add
	gphys = gphys.rename("tr_tppn").
	  set_att("ape_name","precipitation_flux").
	  set_att("units","kg m-2 s-1")
	GPhys::NetCDF_IO.write(file, gphys)
      end

      # PS
      if  File.exist?("PS.ctl") then
	print "GrADS to NetCDF PS\n"
	@data = Ape.new("#{grfile_path}/#{num}/PS.ctl")
	gphys = @data.gphys_open("PS").cut(true,-30..30,1,true)
	gphys = gphys.rename("tr_mslp").
	  set_att("ape_name","air_pressure_at_sea_level").
	  set_att("units","Pa")
	GPhys::NetCDF_IO.write(file, gphys)
	
      end

      # U
      if  File.exist?("U.ctl") then
	print "GrADS to NetCDF U\n"
	@data = Ape.new("#{grfile_path}/#{num}/U.ctl")

	gphys = @data.gphys_open("U").cut(true,-30..30,0.25,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_u250s").
	  set_att("ape_name","eastward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

	gphys = @data.gphys_open("U").cut(true,-30..30,0.570,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_u570s").
	  set_att("ape_name","eastward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

	gphys = @data.gphys_open("U").cut(true,-30..30,0.85,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_u850s").
	  set_att("ape_name","eastward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

	gphys = @data.gphys_open("U").cut(true,-30..30,0.95,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_u950s").
	  set_att("ape_name","eastward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

      end

      # V
      if  File.exist?("V.ctl") then
	print "GrADS to NetCDF V\n"
	@data = Ape.new("#{grfile_path}/#{num}/V.ctl")

	gphys = @data.gphys_open("V").cut(true,-30..30,0.25,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_v250s").
	  set_att("ape_name","northward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

	gphys = @data.gphys_open("V").cut(true,-30..30,0.570,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_v570s").
	  set_att("ape_name","northward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

	gphys = @data.gphys_open("V").cut(true,-30..30,0.85,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_v850s").
	  set_att("ape_name","northward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

	gphys = @data.gphys_open("V").cut(true,-30..30,0.95,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_v950s").
	  set_att("ape_name","northward_wind").
	  set_att("units","m s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

      end

      # Z
      if  File.exist?("Z.ctl") then
	print "GrADS to NetCDF Z\n"
	@data = Ape.new("#{grfile_path}/#{num}/Z.ctl")

	gphys = @data.gphys_open("Z").cut(true,-30..30,0.25,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_phi250s").
	  set_att("ape_name","geopotential_height").
	  set_att("units","m"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

	gphys = @data.gphys_open("Z").cut(true,-30..30,0.85,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_phi850s").
	  set_att("ape_name","geopotential_height").
	  set_att("units","m"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)
      end

      # T
      if  File.exist?("T.ctl") then
	print "GrADS to NetCDF T\n"
	@data = Ape.new("#{grfile_path}/#{num}/T.ctl")
	gphys = @data.gphys_open("T").cut(true,-30..30,0.57,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_t570s").
	  set_att("ape_name","air_temperature").
	  set_att("units","K"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

      end

      # Q
      if  File.exist?("Q.ctl") then
	print "GrADS to NetCDF Q\n"
	@data = Ape.new("#{grfile_path}/#{num}/Q.ctl")

	gphys = @data.gphys_open("Q").cut(true,-30..30,0.85,true)
	lost_axes = gphys.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")
	gphys = gphys.rename("tr_q850s").
	  set_att("ape_name","specific humidity").
	  set_att("units","1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)

      end

      file.close
    }

      Dir.chdir(pwd_path)

  end

end