#!/usr/bin/env ruby
=begin
= mknc_25means.rb -- 年/季節平均 気候値 netCDF ファイル.

以下の量の帯状平均量を保存.

   (a)オイラー平均質量流線関数,
   (b)東西風,
   (c)温度,
   (d)比湿,
   (e)放射収支(大気上端正味短波フラックス,
               大気上端上向き長波フラックス,
               地表面正味短波フラックス,
               地表面正味長波フラックス),
   (f)熱収支(地表面正味潜熱フラックス,
             地表面正味顕熱フラックス,
             降水),
   (g)南北乾燥静的エネルギー輸送量
      (合計, 帯状平均成分, 帯状擾乱成分),
   (h)南北潜熱輸送量
      (合計, 帯状平均成分, 帯状擾乱成分),
   (i)EP フラックスとその発散
      (EP フラックス, 発散),
   (j)残差循環の質量流線関数 を一枚にまとめた図.

== オプション

  1979 - 2003 年 の年平均値
    $ ruby mknc_25means.rb --output hoge.nc --kikan annual --range 1979:2003

  1979 年 の冬平均値
    $ ruby mknc_25means.rb --kikan winter --range 1979:1979
  ...
  

=end

require "getopts"        # for option_parse
require "numru/ggraph"
require "libgphys-n"
include NumRu
include Misc::EMath


##-----------------------
## 25 年平均のオブジェクトを作成

def make_mean_gphys(phys, src, season, year, var=phys)
  pathary = []
  y1 = year[0]; y2 = year[1]
#  (M. Ishiwatari,  2005/03/22)  path 変更
#  path = "../#{phys}.#{src}/"
#  path = "../../#{phys}.#{src}/"
  path = "/GFD_Dennou_Work3/dcchart/atmos_global/NCEP/#{phys}.#{src}/"
  season.each do |m|
    y1.upto(y2) do |y|
      y = y-1 if m == "12"
      fn = path + "#{phys}.#{y}.#{src}/#{phys}_#{y}-#{m}_#{src}.nc"
      pathary << fn
    end
  end
#  gp     = mean_gphys_date_weight(pathary, phys.downcase)
  gp     = mean_gphys(pathary, var.downcase)
  return gp
end


##################################################################

unless getopts(
               "", 
               "kikan:", # "annual", "winter", "spring", "summer", "fall" を指定
               "range:", # 平均する年を指定. とある一年の平均の場合
               "output:mknc_25means.nc" # 出力するファイル.
               )
  print "#{$0}:illegal options. please exec this program with -h/--help. \n"
  exit 1
end


# make gphys objects

## set kikan
year  =   $OPT_range.split(/:/).collect{|y| y.to_i }

months =  ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]
winter = ["12", "01", "02"];  spring = ["03", "04", "05"]
summer = ["06", "07", "08"];  fall = ["09", "10", "11"]
annual = winter + spring + summer + fall


case $OPT_kikan
when "annual"
  kikan = annual
when "winter"
  kikan = winter
when "spring"
  kikan = spring
when "summer"
  kikan = summer
when "fall"
  kikan = fall
else
  if months.delete($OPT_kikan)
    kikan = [ $OPT_kikan ]
  else
    raise "#{$OPT_kikan} is invalid value."
  end
end

## define Constants
Infty = 1.0/0.0
L = UNumeric.new(2.5008e6, 'J.Kg-1')
Rad = UNumeric.new(6.37E6, 'm')
G =   UNumeric.new(9.81  , 'm.s-2')


## 質量流線関数
  strm = make_mean_gphys('STRM', 'NCEP', kikan, year).cut( 'level'=>1000..100 )
## 東西風速
  uwnd = make_mean_gphys('UWND', 'NCEP', kikan, year).mean('lon').cut( 'level'=>1000..100 )
## 温度場
  temp = make_mean_gphys('T',    'NCEP', kikan, year).mean('lon').cut( 'level'=>1000..100 )
## 比湿
  shum = make_mean_gphys('SHUM', 'NCEP', kikan, year).mean('lon').cut( 'level'=>1000..100 )

## 放射

  ## 放射
  dswrf = make_mean_gphys('DSWRF', 'NCEP', kikan, year).mean('lon')
  uswrf = make_mean_gphys('USWRF', 'NCEP', kikan, year).mean('lon')
  nswrf = uswrf - dswrf                                   # 正味の短波放射フラックス
  nswrf.data.set_att('long_name',  'NSR at TOP')
  ulwrf = make_mean_gphys('ULWRF', 'NCEP', kikan, year).mean('lon')
  ulwrf.data.set_att('long_name',  'OLR at TOP')
  nlwrs = make_mean_gphys('NLWRS', 'NCEP', kikan, year).mean('lon')  # 正味の長波放射フラックス
  nlwrs.data.set_att('long_name',  'NLR at SFC')
  nswrs = make_mean_gphys('NSWRS', 'NCEP', kikan, year).mean('lon')  # 地表面正味の短波放射フラックス
  nswrs.data.set_att('long_name',  'NSR at SFC')
  nsw = (-nswrf + nswrs).abs # 大気の吸収量(短波放射)
  nsw.data.set_att('long_name',  'SW HEAT. by ATOM.')
  nlw = (-ulwrf + nlwrs).abs # 大気の吸収量(長波放射)
  nlw.data.set_att('long_name',  'LW COOL. by ATOM.')
  nrw = (nlw   - nsw).abs   # 大気の吸収量(トータル)
  nrw.data.set_att('long_name',  'RW COOL. by ATOM.')
  radiation = [ nswrf.abs, ulwrf, nswrs.abs, nlwrs, nsw, nlw, nrw ] 

  #### 潜熱フラックス
  lhtfl = make_mean_gphys('LHTFL', 'NCEP', kikan, year).mean('lon')
  lhtfl.data.set_att('long_name', 'net latent heat flux at surface')
  #### 顕熱フラックス
  shtfl = make_mean_gphys('SHTFL', 'NCEP', kikan, year).mean('lon')
  shtfl.data.set_att('long_name', 'net sensible heat flux at surface')
  #### 降水(潜熱に換算したもの)
  prate = make_mean_gphys('PRATE', 'NCEP', kikan, year).mean('lon')
  lheat = prate*L
  lheat.data.set_att('long_name', 'latent heat with precipitation')
  heat = [ lhtfl, shtfl, lheat, nrw ]
  
  ## 乾燥静的エネルギー
  dsefl  =      make_mean_gphys('DSEFL', 'NCEP', kikan, year, 'dsefl')
  dsefl_bar   = make_mean_gphys('DSEFL', 'NCEP', kikan, year, 'dsefl_bar')
  dsefl_dash  = make_mean_gphys('DSEFL', 'NCEP', kikan, year, 'dsefl_dash')
  
  ## 潜熱輸送量
  lhefl  =      make_mean_gphys('LHEFL', 'NCEP', kikan, year, 'lhefl')
  lhefl_bar   = make_mean_gphys('LHEFL', 'NCEP', kikan, year, 'lhefl_bar')
  lhefl_dash  = make_mean_gphys('LHEFL', 'NCEP', kikan, year, 'lhefl_dash')
  
  # get cos
  lat     = lhefl.coord('lat')
  grid    = Grid.new(lhefl.axis('lat'))
  cos_lat = GPhys.new(grid, cos(lat*2*PI/360.0))
  
  
  dse =      dsefl.integrate('level') /G * Rad * cos_lat * 2 * PI * 100
  dse_bar =  dsefl_bar.integrate('level') /G * Rad * cos_lat * 2 * PI * 100
  dse_dash = dsefl_dash.integrate('level') /G * Rad * cos_lat * 2 * PI * 100
  dse.data.set_att('units', 'W')
  dse_bar.data.set_att('units', 'W')
  dse_dash.data.set_att('units', 'W')  
  
  lhe =      lhefl.integrate('level') /G * Rad * cos_lat * 2 * PI * 100
  lhe_bar =  lhefl_bar.integrate('level') /G * Rad * cos_lat * 2 * PI * 100
  lhe_dash = lhefl_dash.integrate('level') /G * Rad * cos_lat * 2 * PI * 100
  lhe.data.set_att('units', 'W')
  lhe.data.set_att('long_name', 'latent heat energy transport (total)')
  lhe_bar.data.set_att('units', 'W')
  lhe_dash.data.set_att('units', 'W')
  
  # 湿潤静的エネルギー
  mse = lhe + dse
  mse.data.set_att('long_name', 'moist static energy transport (total)')

  dry_ene = [ dse, dse_bar, dse_dash, mse ]
  lat_ene = [ lhe, lhe_bar, lhe_dash, mse ]
  
  ## EP フラックス
  epflx_phi = make_mean_gphys('EPFLX', 'NCEP', kikan, year, 'epflx_phi')
  epflx_p   = make_mean_gphys('EPFLX', 'NCEP', kikan, year, 'epflx_p')
  epflx_div = make_mean_gphys('EPFLX', 'NCEP', kikan, year, 'epflx_div')
  ## 残差循環
  strm_rmean = make_mean_gphys('STRM_RMEAN', 'NCEP', kikan, year)
    

####### delete valid_range

[shum, strm, uwnd, temp, 
 dse, dse_bar, dse_dash, 
 lhe, lhe_bar, lhe_dash,
 epflx_phi, epflx_p, epflx_div,
 strm_rmean,
 nswrf,
 ulwrf,
 nswrs,
 nlwrs,
 lhtfl,
 shtfl,
 lheat
].each{|gp|
  gp.data.att_names.each{|nm|
    gp.del_att(nm) if /^valid_/ =~ nm
    gp.del_att(nm) if /^actual_/ =~ nm
  }
}      
            
#############################################################################

## netCDF ファイルに保存

nc = NetCDF::create('300hPa_'+$OPT_output)
GPhys::IO.write( nc, shum )  	 # 比湿
nc.close

nc = NetCDF::create('100hPa_'+$OPT_output)
GPhys::IO.write( nc, strm )  	 # 質量流線関数
GPhys::IO.write( nc, uwnd )  	 # 東西風速
GPhys::IO.write( nc, temp )  	 # 温度場

GPhys::IO.write( nc, dse )       # 乾燥静的エネルギー輸送量(合計)
GPhys::IO.write( nc, dse_bar )   #                         (平均)
GPhys::IO.write( nc, dse_dash )  #                         (擾乱)
GPhys::IO.write( nc, lhe ) 	 # 潜熱輸送量              (合計)
GPhys::IO.write( nc, lhe_bar )   #                         (平均)
GPhys::IO.write( nc, lhe_dash )  #                         (擾乱)

GPhys::IO.write( nc, epflx_phi ) # EP フラックス 緯度成分
GPhys::IO.write( nc, epflx_p   ) #               鉛直(圧力)成分
GPhys::IO.write( nc, epflx_div ) # EP フラックス 発散
GPhys::IO.write( nc, strm_rmean) # 残差子午面循環 質量流線関数
nc.close

nc = NetCDF::create('GAUSSIAN_'+$OPT_output)
GPhys::IO.write( nc, nswrf ) 	 # 大気上端正味の短波放射フラックス 
GPhys::IO.write( nc, ulwrf ) 	 # 大気上端上向き長波放射フラックス 
GPhys::IO.write( nc, nswrs ) 	 # 地表面正味の短波放射フラックス
GPhys::IO.write( nc, nlwrs ) 	 # 地表面正味の長波放射フラックス
GPhys::IO.write( nc, lhtfl ) 	 # 潜熱フラックス
GPhys::IO.write( nc, shtfl ) 	 # 顕熱フラックス
GPhys::IO.write( nc, lheat ) 	 # 降水
nc.close