#!/usr/bin/env ruby
=begin
= mkfig_25time_series.rb -- make time-series figure of zonal wind intensity

=end

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

def cal_correlation(x, y) # 相関係数を求める
  xa = x - x.mean(0)           # x の偏差
  ya = y - y.mean(0)           # y の偏差
  xstddev = x.stddev           # x の標準偏差(分散の平方根)
  ystddev = y.stddev           # y   ::
  
  crc = (( xa*ya ).mean(0))/(xstddev*ystddev) # 相関係数 
  return crc
end



################################################################
#                        make gphys 
################################################################
year = [1979, 2003]
#year = [1979, 1980]
month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]   # 初期値 
#month = ["02", "08"]   # 初期値 

uwnd_nh_20_40_long_name = "mean zonal wind speed (N20 - N40)"
uwnd_nh_40_60_long_name = "mean zonal wind speed (N40 - N60)"
uwnd_sh_20_40_long_name = "mean zonal wind speed (S20 - S40)"
uwnd_sh_40_60_long_name = "mean zonal wind speed (S40 - S60)"

## make compisit-GPhys object 
# make array which set the gphys objects.

axmonth = []
i = 0

uwnd_name = ""
uwnd_unit = ""

uwnd_north_20_40 = []
uwnd_north_40_60 = []
uwnd_south_20_40 = []
uwnd_south_40_60 = []

year[0].upto(year[1]) do |y|
  month.each do |m|
    path = "../../../UWND.NCEP/UWND.#{y}.NCEP/UWND_#{y}-#{m}_NCEP.nc" # stream function
    gp_uwnd = GPhys::NetCDF_IO.open(path, 'uwnd').cut('level'=>100..1000).mean('lon')
    gp_climat = GPhys::NetCDF_IO.open("../../25MEANS.NCEP/25MEANS.1979-2003.NCEP/100hPa_PHYS_1979-2003_#{m}_NCEP.nc", 'uwnd')
#    if m == "02"
      uwnd_north_20_40 <<  (gp_uwnd).cut('level'=>200,  'lat'=>20..40  ).mean('lat').val - (gp_climat).cut('level'=>200,  'lat'=>20..40  ).mean('lat').val
      uwnd_north_40_60 <<  (gp_uwnd).cut('level'=>200,  'lat'=>40..60  ).mean('lat').val - (gp_climat).cut('level'=>200,  'lat'=>40..60  ).mean('lat').val
#    elsif m == "08"
      uwnd_south_20_40 <<  (gp_uwnd).cut('level'=>200,  'lat'=>-20..-40).mean('lat').val - (gp_climat).cut('level'=>200,  'lat'=>-20..-40).mean('lat').val
      uwnd_south_40_60 <<  (gp_uwnd).cut('level'=>200,  'lat'=>-40..-60).mean('lat').val - (gp_climat).cut('level'=>200,  'lat'=>-40..-60).mean('lat').val
#    end
    uwnd_name =      gp_uwnd.data.name
    uwnd_unit =      gp_uwnd.data.get_att("units")
  end
end


### for intensity

  na_uwnd_north_20_40 = NArray.to_na(uwnd_north_20_40)
  na_uwnd_north_40_60 = NArray.to_na(uwnd_north_40_60)
  na_uwnd_south_20_40 = NArray.to_na(uwnd_south_20_40)
  na_uwnd_south_40_60 = NArray.to_na(uwnd_south_40_60)


############      stream function        #########################

strm_nh_hadley_long_name = "hadley circulation intencity (NH)"
strm_sh_hadley_long_name = "hadley circulation intencity (SH)"

## make compisit-GPhys object 
# make array which set the gphys objects.

axmonth = []
i = 0

strm_name = ""
strm_unit = ""

strm_hadley_north = []
strm_hadley_south = []

year[0].upto(year[1]) do |y|
  month.each do |m|
    path = "../../../STRM.NCEP/STRM.#{y}.NCEP/STRM_#{y}-#{m}_NCEP.nc" # stream function
    gp_climat = GPhys::NetCDF_IO.open("../../25MEANS.NCEP/25MEANS.1979-2003.NCEP/100hPa_PHYS_1979-2003_#{m}_NCEP.nc", 'strm')
    gp_strm =   GPhys::NetCDF_IO.open(path, 'strm').cut('level'=>100..1000) 

    gp_strm = GPhys::NetCDF_IO.open(path, 'strm')
#    if m == "02"
      strm_hadley_north <<  gp_strm.cut( 'lat'=>  0..30  ).data.val.max - gp_climat.cut( 'lat'=>  0..30  ).data.val.max
#    elsif m == "08"
      strm_hadley_south <<  gp_strm.cut( 'lat'=>  0..-30 ).data.val.min - gp_climat.cut( 'lat'=>  0..-30 ).data.val.min
#    end
    strm_name =      gp_strm.data.name
    strm_unit =      gp_strm.data.get_att("units")

  end
end


### for intensity

  na_strm_hadley_north = NArray.to_na(strm_hadley_north)
  na_strm_hadley_south = NArray.to_na(strm_hadley_south)

                   
################################################################
#                        描画ルーチン
################################################################



##
# 事前準備

DCL.uscset('cyspos', 'B' )              # y 軸の単位の位置を下方へ 
rsizel2 = DCL.uzrget('rsizel2')         # 現在のラベルサイズを取得
DCL.uzrset('rsizel2', rsizel2*0.42 )    # ラベルサイズをデフォルトの 0.5 倍に

##
# お絵描きメイン

DCL.swpset('lsep',  true)    # ページ毎別々のファイルに落す
DCL.gropn(2)
DCL.sldiv('T', 2, 2)
DCL.sgpset('lcntl', false)   # 
DCL.sgpset('lfull',true)     # フルスクリーン
DCL.uzfact(0.55)             # 
DCL.sgpset('lcorner',false)  # コーナーを取っ払う 
DCL.sgpset('lfprop',true)    # 
DCL.udpset('lmsg',false)     # 
DCL.uscset('cyspos', 'B' )   # y 軸の単位の位置を下方へ 

## y 軸範囲

axx = NArray[-5e10, 5e10] # y 軸範囲. -7m/s から + 7m/s
axy = NArray[-7, 7]       # y 軸範囲. -7m/s から + 7m/s


### 1 枚目 : 北半球 (N20 - N40)

DCL::grfrm                   # ページ確定 

x = na_strm_hadley_north
y = na_uwnd_north_20_40

DCL::usspnt(axx, axy)
DCL::uspfit
DCL::grstrf

DCL::ussttl('Deviation of mass stream function (NH)', strm_unit, 'Deviation of zonal wind ', uwnd_unit )
DCL::usdaxs
DCL::uxmttl('T', 'Zonal Wind at 200hPa (N20 - N40)', 0.0)

DCL::uusmkt(2)
DCL::uumrk(x, y)

print "北半球ハドレー循環強度と200 hPa 面東西風(N20-N40)", cal_correlation(x, y), "\n"


### 2 枚目 : 北半球 (N40 - N60)

DCL::grfrm                   # ページ確定 

x = na_strm_hadley_north
y = na_uwnd_north_40_60

DCL::usspnt(axx, axy)
DCL::uspfit
DCL::grstrf

DCL::ussttl('Deviation of mass stream function (NH)', strm_unit, 'Deviation of zonal wind ', uwnd_unit )
DCL::usdaxs
DCL::uxmttl('T', 'Zonal Wind at 200hPa (N40 - N60)', 0.0)

DCL::uusmkt(2)
DCL::uumrk(x, y)

print "北半球ハドレー循環強度と200 hPa 面東西風(N40-N60)", cal_correlation(x, y), "\n"

### 3 枚目 : 南半球 (S20 - S40)

DCL::grfrm                   # ページ確定 

x = na_strm_hadley_south
y = na_uwnd_south_20_40

DCL::usspnt(axx, axy)
DCL::uspfit
DCL::grstrf

DCL::ussttl('Deviation of mass stream function (SH)', strm_unit, 'Deviation of zonal wind ', uwnd_unit )
DCL::usdaxs
DCL::uxmttl('T', 'Zonal Wind at 200hPa (S20 - S40)', 0.0)

DCL::uusmkt(2)
DCL::uumrk(x, y)

print "南半球ハドレー循環強度と200 hPa 面東西風(S20-S40)", cal_correlation(x, y), "\n"

### 4 枚目 : 南半球 (S40 - S60)

DCL::grfrm                   # ページ確定 

x = na_strm_hadley_south
y = na_uwnd_south_40_60

DCL::usspnt(axx, axy)
DCL::uspfit
DCL::grstrf

DCL::ussttl('Deviation of mass stream function (SH)', strm_unit, 'Deviation of zonal wind ', uwnd_unit )
DCL::usdaxs
DCL::uxmttl('T', 'Zonal Wind at 200hPa (S40 - S60)', 0.0)

DCL::uusmkt(2)
DCL::uumrk(x, y)

print "南半球ハドレー循環強度と200 hPa 面東西風(S40-S60)", cal_correlation(x, y), "\n"

=begin
# マーカー情報記述
charsize = 1.0 * DCL.uzpget('rsizec1')
dvx = 0.01
dvy = charsize*1.5
vxmin,vxmax,vymin,vymax = DCL.sgqvpt
vx = 0.12
vy = 0.1 + charsize/2
DCL::sgpmu([vx], [vy])
DCL::sgtxzv(vx+dvx,vy,uwnd_nh_20_40_long_name, 
	    charsize, 0, -1, 1)
=end
DCL.grcls

##################################################
# 画像ファイル名を変更
File.rename('dcl_001.ps', 'UWND_1973-2003_TIME-SERIES_NCEP_DEVIATION_area_mean_sanpuzu.ps')