#!/usr/bin/env ruby

=begin

=gpcat

複数ファイルに対して時空間方向に接続して X 平均を取る. 
指定されたディレクトリ以下の time* というディレクトリに対して再帰的に実行. 

==USAGE

  % arare-addmean-velz.rb --dir Directory


==OPTIONS

      --help         : print this message. 
      --dir dirvar   : directory name (required).

==HISTORY

  2009/04/23  K Sugiyama (created)

=end

require "numru/netcdf"
require "numru/dcl"
require "numru/ggraph"
require "getoptlong"
include NumRu

###
###引数処理
###
parser = GetoptLong.new
parser.set_options(
                   ###    global option   ###
                   ['--dir',  "-d", GetoptLong::REQUIRED_ARGUMENT],
                   ["--help", "-h", GetoptLong::NO_ARGUMENT ]
                   )
begin
  parser.each_option do |name, arg|
    eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'" 
  end
    rescue
  exit(1)
end

if $OPT_HELP then
  help
  exit(0)
end  

###
### 初期設定
###
coords = ["x", "y", "z", "t", "s", "alt", "t_nv", "t_bnds"] # 軸の設定
dir = "#{$OPT_dir}"           # トップディレクトリ

# ファイル置き場
meanfiles = "meanfiles"	 
submeanfiles = "#{meanfiles}"
Dir::mkdir(meanfiles) unless FileTest.exist?(meanfiles)
Dir::mkdir(submeanfiles) unless FileTest.exist?(submeanfiles)

###
### 水平平均を行う
###   水平平均された量に関しては, ファイル名を変更するだけ.
###

Dir.glob("#{dir}/time*").each{ |subdir|

  output1 = "#{submeanfiles}/#{File::basename(subdir)}_VelZRMS.nc"
  output2 = "#{submeanfiles}/#{File::basename(subdir)}_VelZSKW.nc"
  output3 = "#{submeanfiles}/#{File::basename(subdir)}_VelZKRT.nc"

  Dir.glob("#{subdir}/*all*VelZ.nc").each{ |ncfile|
    gturl = "#{ncfile}@VelZ"
    p gturl
    gphys = GPhys::IO.open_gturl(gturl)
      
    MEAN = gphys.mean('x')
    DV   = gphys - MEAN
    RMS  = ((gphys ** 2.0e0).mean('x').mean('y')).sqrt
    SDV  = ((DV ** 2.0e0).mean('x').mean('y')).sqrt
    SKW  = ((DV ** 3.0e0).mean('x').mean('y')) / (SDV ** 3.0e0)
    KRT  = ((DV ** 4.0e0).mean('x').mean('y')) / (SDV ** 4.0e0) - 3.0e0
      
    outncfile=NetCDF.create(output1)
    GPhys::IO.write( outncfile, RMS.rename('VelZRMS') ) # Output GPhys variable
    outncfile.close
    p "#{gturl} is written to #{output1}."

    outncfile=NetCDF.create(output2)
    GPhys::IO.write( outncfile, SKW.rename('VelZSKW') ) # Output GPhys variable
    outncfile.close
    p "#{gturl} is written to #{output2}."
    
    outncfile=NetCDF.create(output3)
    GPhys::IO.write( outncfile, KRT.rename('VelZKRT') )  # Output GPhys variable
    outncfile.close
    p "#{gturl} is written to #{output3}."
  }
}

###
### 時間方向に接続する. 
### ファイル名が正しく時系列に並んでいることを前提とする. 
###
vars = ["VelZRMS", "VelZSKW", "VelZKRT"]
prefix = "#{submeanfiles}/AllData"

vars.each{ |var|
  output = "#{prefix}_#{var}.nc"
  ncfiles2 = Dir::glob("#{submeanfiles}/*_#{var}.nc").sort
  p ncfiles2, output, var
  unless ncfiles2.size == 1 
    gphys3 = GPhys::IO.open(ncfiles2, var)
    print "GPhys variable '#{var}' in NetCDF files, " + ncfiles2.join(', ') +", was opened successfully.\n"
  else
    gphys3 = GPhys::IO.open(ncfiles2[0], var)
    print "GPhys variable '#{var}' in NetCDF files, " + ncfiles2[0] +", was opened successfully.\n"
  end
  
  outncfile=NetCDF.create(output)
  GPhys::IO.write( outncfile, gphys3 )
  p "#{var} is written to #{output}."
  outncfile.close
} 


