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

# $rezol = "T159L48_eml"
# 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,true,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/"
    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 = ["0011"]
    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

#      # OLR
#      if  File.exist?("OLR.ctl") then
#	print "GrADS to NetCDF OLR\n"
#	@data = Ape.new("#{grfile_path}/#{num}/OLR.ctl")
#	gphys = @data.gphys_open("OLR")[true,true,0,true]
#	gphys = gphys.rename("tr_lw_toa").
#	  set_att("ape_name","toa_outgoing_longwave_flux").
#	  set_att("units","W m-2")
#	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,true,25000,true)
#	lost_axes = gphys.lost_axes.to_s.sub("z=","plev=")
#	gphys = gphys.rename("Z").
#	  set_att("ape_name","geopotential_height").
#	  set_att("units","m"). 
#	  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","geopotential_height").
	  set_att("units","m"). 
	  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","1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)
      end


      # 仮温度
      # Tv=T(1+0.61q)
      if  File.exist?("Q.ctl") then
	print "GrADS to NetCDF Tv \n"
	@data = Ape.new("#{grfile_path}/#{num}/T.ctl")
	gphys_t = @data.gphys_open("T").cut(true,0,true,true)
	@data = Ape.new("#{grfile_path}/#{num}/Q.ctl")
	gphys_q = @data.gphys_open("Q").cut(true,0,true,true)
	lost_axes = gphys_t.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys_t*(1+0.61*gphys_q.data.val)
	gphys = gphys.rename("TV").
	  set_att("ape_name","virtual temperature").
	  set_att("units","K"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)
      end


      # pf_t_conv
      if  File.exist?("DTCUM.ctl") && File.exist?("DTLSC.ctl") then
	print "GrADS to NetCDF DTCUM & DTLSC\n"
	@data = Ape.new("#{grfile_path}/#{num}/DTCUM.ctl")
	gphys = @data.gphys_open("DTCUM").cut(true,0,true,true)
	gphys_add = Ape.new("#{grfile_path}/#{num}/DTLSC.ctl").
	  gphys_open("DTLSC").cut(true,0,true,true)
	gphys = gphys + gphys_add
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("pf_t_conv").
	  set_att("ape_name",
		  "tendency_of_air_temperature_due_to_moist_convection").
	  set_att("units","K s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)
      end

      # DTRADL
      if  File.exist?("DTRADL.ctl") then
	print "GrADS to NetCDF DTRADL\n"
	@data = Ape.new("#{grfile_path}/#{num}/DTRADL.ctl")
	gphys = @data.gphys_open("DTRADL").cut(true,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("DTRADL").
	  set_att("ape_name",
		  "tendency_of_air_temperature_due_to_longwave_heating").
	  set_att("units","K s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)
      end

      # DTRADS
      if  File.exist?("DTRADS.ctl") then
	print "GrADS to NetCDF DTRADS\n"
	@data = Ape.new("#{grfile_path}/#{num}/DTRADS.ctl")
	gphys = @data.gphys_open("DTRADS").cut(true,0,true,true)
	lost_axes = gphys.lost_axes.to_s.sub("y=","lat=")
	gphys = gphys.rename("DTRADS").
	  set_att("ape_name",
		  "tendency_of_air_temperature_due_to_shortwave_heating").
	  set_att("units","K s-1"). 
	  set_att("lost_axis",lost_axes)
	GPhys::NetCDF_IO.write(file, gphys)
      end

      file.close
    }

      Dir.chdir(pwd_path)

  end

end