#!/usr/bin/env ruby

=begin

表題: AGUforAPE GrADS to NetCDF

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

=end

END{

mknetcdf

}

# ----------------------------------------------
# 実験設定

$sstid = "expID#{ARGV[0]}"
$cumulus = "non"
$resol = "T79L48_#{$cumulus}"
$expdir = "/work11/ape/grads/#{$resol}_#{$sstid}"
$outdir = "/work11/ape/grads/mean"
#$outdir = "/work11/ape/yukiko"

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

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

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

$local_path = '/work11/ape/yukiko/lib'
$: << $local_path

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

require "#{$local_path}/ape-view.rb"

# -----------------------------------------------
# Ape 変数リスト

$cntfile_array = [
  "U",
  "V",
  "T",
  "OMG",
  "Z",
  "Q",
  "QL",
  "RH",
  'PS',
  'OSRD',
  "OSRU",
  "OLR",
  "CCOVER",
  "PRCPC",
  "PRCPL",
  "SSRD",
  "SSRU",
  "SLRD",
  "SLRU",
  "EVAP",
  "SENS",
  "TAUX",
  "TAUY",
  "VICLDW",
  "CLDW",
  "DTDYN",
  "DTRADS",
  "DTRADL",
  "DTVDF",
  "DTCUM",
  "DTLSC",
  "DTDCA",
  "DQDYN",
  "DQVDF",
  "DQLSC",
  "DQDCA",
  "SRMNQ",
  "DUDYN",
  "DUVDF",
  "DUGWD",
  "DVDYN",
  "DVVDF",
  "DVGWD",
  "SLP"
  ]

$unit_hash = {
  'U'      =>'m/s',
  'V'      =>'m/s',
  'T'      =>'K',
  'OMG'    =>'Pa/s',
  'Z'      =>'m',
  'Q'      =>'kg/kg',
  'QL'     =>'kg/kg',
  'RH'     =>'1',
  'PS'     =>'Pa',
  'OSRD'   =>'W/m**2',
  'OSRU'   =>'W/m**2',
  'OLR'    =>'W/m**2',
  'CCOVER' =>'1',
  'PRCPC'  =>'kg/m**2/s',
  'PRCPL'  =>'kg/m**2/s',
  'SSRD'   =>'W/m**2',
  'SSRU'   =>'W/m**2',
  'SLRD'   =>'W/m**2',
  'SLRU'   =>'W/m**2',
  'EVAP'   =>'W/m**2',
  'SENS'   =>'W/m**2',
  'TAUX'   =>'N/m**2',
  'TAUY'   =>'N/m**2',
  'VICLDW' =>'kg/m**2',
  'CLDW'   =>'kg/kg',
  'DTDYN'  =>'K/s',
  'DTRADS' =>'K/s',
  'DTRADL' =>'K/s',
  'DTVDF'  =>'K/s',
  'DTCUM'  =>'K/s',
  'DTLSC'  =>'K/s',
  'DTDCA'  =>'K/s',
  'DQDYN'  =>'1/s',
  'DQVDF'  =>'1/s',
  'DQCUM'  =>'1/s',
  'DQLSC'  =>'1/s',
  'DQDCA'  =>'1/s',
  'SRMNQ'  =>'1/s',
  'DUDYN'  =>'m/s**2',
  'DUVDF'  =>'m/s**2',
  'DUGWD'  =>'m/s**2',
  'DVDYN'  =>'m/s**2',
  'DVVDF'  =>'m/s**2',
  'DVGWD'  =>'m/s**2',
  'SLP'   =>'Pa'
}

$longname_hash = {
  'U'      =>'u-velocity',
  'V'      =>'v-velocity',
  'T'      =>'temperature',
  'OMG'    =>'p-velocity',
  'Z'      =>'geopotential height',
  'Q'      =>'specific humidity',
  'QL'     =>'liquid water',
  'RH'     =>'relative humidity',
  'PS'     =>'surface pressure',
  'OSRD'   =>'top.short.down',
  'OSRU'   =>'top.short.up',
  'OLR'    =>'top longwave',
  'CCOVER' =>'cloud cover',
  'PRCPC'  =>'cumulus precipitation',
  'PRCPL'  =>'L.S.precipitation',
  'SSRD'   =>'surf.short.down',
  'SSRU'   =>'surf.short.up',
  'SLRD'   =>'surf.long.down',
  'SLRU'   =>'surf.long.up',
  'EVAP'   =>'latent heat flux',
  'SENS'   =>'sensible heat flux',
  'TAUX'   =>'surface stress(x)',
  'TAUY'   =>'surface stress(y)',
  'VICLDW' =>'integrated cldw',
  'CLDW'   =>'cloud water',
  'DTDYN'  =>'dynamics T-tendency',
  'DTRADS' =>'radiative heating(short)',
  'DTRADL' =>'radiative heating(long)',
  'DTVDF'  =>'diffusion heating',
  'DTCUM'  =>'cumulus heating',
  'DTLSC'  =>'L.S.cond. heating',
  'DTDCA'  =>'dry conv. heating',
  'DQDYN'  =>'dynamics q-tendency',
  'DQVDF'  =>'diffusion moistning',
  'DQCUM'  =>'cumulus moistning',
  'DQLSC'  =>'L.S.cond. moistening',
  'DQDCA'  =>'dry conv. moistning',
  'SRMNQ'  =>'moisture fixer',
  'DUDYN'  =>'dynamics u-tendency',
  'DUVDF'  =>'diffusion du/dt',
  'DUGWD'  =>'gravity wave du/dt',
  'DVDYN'  =>'dynamics v-tendency',
  'DVVDF'  =>'diffusion dv/dt',
  'DVGWD'  =>'gravity wave dv/dt',
  'SLP'   =>'sea level pressure'
}

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

def mknetcdf
  
  # ファイル一覧取得
#  $cntfile_array = []
  Dir.foreach("#{$expdir}/mean") { |item|
    if item =~ /ctl$|CTL$/
      puts item
#      $cntfile_array.push( File.basename(item, ".ctl").sub("mean","") ) 
    end
  }
  
  # netCDF 初期化
  file = NetCDF.create("#{$outdir}/#{$resol}_#{$sstid}_mean.nc")

  # 書き込み
  $cntfile_array.each {|item|

    if  File.exist?("#{$expdir}/mean/mean#{item}.ctl") then
      data = Ape.new("#{$expdir}/mean/mean#{item}.ctl")
    
      if item =~ /CCOVER$|EVAP$|OLR$|OSRD$|OSRU$|PRCPC$|PRCPL$|PS$|SENS$|SLP$|SLRD$|SLRU$|SSRD$|SSRU$|TAUX$|TAUY$|VICLDW$/
	print "#{item} -- GrADS to NetCDF -- #{Time.now} ** lev1 \n"
	if $resol =~/T159/
	  gphys = data.gphys_open(item)[true,true,0,0]
	else
	  gphys = data.gphys_open(item)[true,true,0,true].mean(2)
	end
      else
	print "#{item} -- GrADS to NetCDF -- #{Time.now} \n"
	if $resol =~/T159/
	  gphys = data.gphys_open(item)[true,true,true,0]
	else
	  gphys = data.gphys_open(item).mean(3)
	end
      end

      gphys = gphys. 
	set_att("long_name",$longname_hash[item]).
	set_att("units",$unit_hash[item]). 
	dimname_change
      GPhys::NetCDF_IO.write(file, gphys)
   end
  }
  
  file.put_att("Conventions", "CF-1.0")
  file.put_att("title","Aqua Planet: #{$resol}_#{$sstid} Experiment")
  file.put_att("history", "Original data produced: #{Time.now} from 3 years mean GrADS files on ES")
  file.put_att("institution", "AGU for APE")
  file.put_att("source","afes_115_ape_#{$cumulus}.LM")

  file.close
  
end


class GPhys

  def dimname_change

    xdim = VArray.new(@grid.axis(0).pos.val).rename("lon").
      put_att("standard_name","longitude").
      put_att("units","degrees_east")
    ydim = VArray.new(@grid.axis(1).pos.val).rename("lat"). 
      put_att("standard_name","latitude").
      put_att("units","degrees_north")

    if @data.val.dim  == 3
      zdim = VArray.new(@grid.axis(2).pos.val).rename("sigma").
	put_att("standard_name","sigma_level").
	put_att("units","1")
    end

    xdim = Axis.new().set_pos(xdim)
    ydim = Axis.new().set_pos(ydim)
    zdim = Axis.new().set_pos(zdim)       if @data.val.dim == 3

    @grid = @grid.change_axis(0, xdim)
    @grid = @grid.change_axis(1, ydim)
    @grid = @grid.change_axis(2, zdim)    if @data.val.dim == 3
 
#    @grid = Grid.new(xdim,ydim,zdim)
 
    # gphys class の値を返す
    GPhys.new( @grid, @data )

  end

end