#!/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 = "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

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

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

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

 $rezol = "T159L48_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


}


# ----------------------------------------------
# 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 mkprcpt
  
    grdir = Array.new
    Dir.foreach($grfile_path) { |item|
      if item =~ /^00/ 
	grdir.push(item)
	print "#{item}\n"
      end
    }
    
    # [PRCPC, PRCPL, PS]
    grdir.each{ |num|
      
      @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
      gphys = gphys.rename("PRCPT").
	set_att("ape_name","precipitation flux").
	set_att("units","kg m-2 s-1")

      file = NetCDF.create("#{$gr2ncfile_path}/#{$rezol}_expID01_#{num}.nc")
      GPhys::NetCDF_IO.write(file, gphys)
      file.close
    }
    
  end

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

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

    pwd_path = Dir.pwd

#    grdir= ["0001","0002"]
#    grdir= ["0003"]
    grdir.each{ |num|

      file = NetCDF.create("#{gr2ncfile_path}/#{$rezol}_expID01_#{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,true,1,true)
	gphys_add = Ape.new("#{grfile_path}/#{num}/PRCPL.ctl").
	  gphys_open("PRCPL").cut(true,true,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,true,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,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("U").
	  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,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("V").
	  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,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("Z").
	  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

      # OMG
      if  File.exist?("OMG.ctl") then
	print "GrADS to NetCDF OMG\n"
	@data = Ape.new("#{grfile_path}/#{num}/OMG.ctl")
	gphys = @data.gphys_open("OMG").cut(true,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("OMG").
	  set_att("ape_name","omega").
	  set_att("units","Pa s-1"). 
	  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,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("T").
	  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,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("Q").
	  set_att("ape_name","specific humidity").
	  set_att("units","kg/kg"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)
      end

      file.close
    }


      Dir.chdir(pwd_path)

  end

end