#! /usr/bin/ruby
# -*- coding: utf-8 -*-
#  乾燥静的エネルギー南北輸送, 潜熱南北輸送量
#


# Variables
#   gp_DryStatEngy : Dry Static Energy
#   gp_LatentEngy  : Latent Energy


# 色つけるかどうかのオプションがあると良い?

# 使用方法
#  % fluxmeri.pl --time_start=366 --time_end=730 --wsn=2


require "getoptlong"        # for option_parse
require "numru/ggraph"
require "pp"
include NumRu


parser = GetoptLong.new
parser.set_options(
  ###    global option   ###
  ['--time_start',                GetoptLong::REQUIRED_ARGUMENT],
  ['--time_end',                  GetoptLong::REQUIRED_ARGUMENT],
  ['--max',                GetoptLong::REQUIRED_ARGUMENT],
  ['--min',                GetoptLong::REQUIRED_ARGUMENT],
  ['--wsn',                       GetoptLong::REQUIRED_ARGUMENT]
)
parser.each_option do |name, arg|
    eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'"  # strage option value to $OPT_val
end

#
# Settings for option values
#
time_start = ($OPT_time_start||1).to_f
time_end = ($OPT_time_end||9999).to_f
min = ($OPT_min||-1.0e12).to_f
max = ($OPT_max||1.0e12).to_f

#
# Constants
#
LatentHeat = 2.5e6
CpDry = 1004.6       # $ C_p $ [J kg-1 K-1]
Grav = 9.8
RPlanet = 6.371 * 10 ** 6  # [m]

# Phi の作成
if File.exist?('Phi.nc') then
  system("rm Phi.nc")
end

p 'generating Phi.nc...'
system("ruby /GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/geopotential_time.rb --var Temp.nc@Temp --range_time #{time_start}:#{time_end}")


# Data Input
# meridional velocity
var = 'V'
path = var + '.nc'
gp_v = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)
# temperature
var = 'Temp'
path = var + '.nc'
gp_temp = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)
# geopotential hight
var = 'Phi'
path = var + '.nc'
gp_phi = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)
# specific humidity
var = 'QVap'
path = var + '.nc'
gp_qvap = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)

# latitude grid
lat = GPhys::IO.open(path, 'lat')
lat_rad = lat * Math::PI/180.0
LatLength = 2 * Math::PI * RPlanet * lat_rad.cos

# sigma weight
sig_weight = GPhys::IO.open(path, 'sig_weight')
na_sig_weight = sig_weight.val

nsig = na_sig_weight.size


#
# Produce necessary data
# 外で作るようにした方が良いかも.
#
gp_DryStatEngy = gp_temp * CpDry + gp_phi
gp_LatentEngy = gp_qvap * LatentHeat


#
# Data processing
#

# gt3 時代は解析期間だけ時間方向に切り出したファイルとして
# $file-anal.out を用意していた.


# _
# v(x,y,z) : Time mean
# gtool 時代は $file-avr.out
gp_v_tmean = gp_v.mean('time')
gp_DryStatEngy_tmean = gp_DryStatEngy.mean('time')
gp_LatentEngy_tmean = gp_LatentEngy.mean('time')

#  _
# [v](y,z) : Time mean, Zonal mean
# gt3 の時代には $file-avr-zm.out
gp_v_tmean_xmean = gp_v_tmean.mean('lon')
gp_DryStatEngy_tmean_xmean = gp_DryStatEngy_tmean.mean('lon')
gp_LatentEngy_tmean_xmean = gp_LatentEngy_tmean.mean('lon')


#                               _
# sigma mean of [v](y,z) : int [v] dσ
# gt3 の時代には $file-avr-zm-zmean.out
gp_v_tmean_xmean_zmean = (gp_v_tmean_xmean * sig_weight).sum('sig')
gp_DryStatEngy_tmean_xmean_zmean = (gp_DryStatEngy_tmean_xmean * sig_weight).sum('sig')
gp_LatentEngy_tmean_xmean_zmean = (gp_LatentEngy_tmean_xmean * sig_weight).sum('sig')

#                              _         _
# eddy comonent of [v](y,z) : [v] - int [v] dσ
# gt3 の時代には $file-avr-zm-zeddy.out
gp_v_tmean_xmean_zeddy = gp_v_tmean_xmean - gp_v_tmean_xmean_zmean
gp_LatentEngy_tmean_xmean_zeddy = gp_LatentEngy_tmean_xmean - gp_LatentEngy_tmean_xmean_zmean
gp_DryStatEngy_tmean_xmean_zeddy = gp_DryStatEngy_tmean_xmean - gp_DryStatEngy_tmean_xmean_zmean





#  _
# 時間平均場の東西平均からのずれ v*(x,y,z) 
# gt3 時代は $file-avr-eddy.out 
gp_v_tmean_xeddy = gp_v_tmean.eddy('lon')
gp_DryStatEngy_tmean_xeddy = gp_DryStatEngy_tmean.eddy('lon')
gp_LatentEngy_tmean_xeddy = gp_LatentEngy_tmean.eddy('lon')


#
# v'(t,z,y,z) (で良いのかな?) : 時間方向からのずれ
# gt3 時代は $file-dist.out
gp_v_tdist = gp_v.eddy('time')
gp_DryStatEngy_tdist = gp_DryStatEngy.eddy('time')
gp_LatentEngy_tdist = gp_LatentEngy.eddy('time')


#
# 
# Calculations for energy transport
#    Loop structure is more better?
#
#


#  __
# [Sv] : Total Transport
#  gt3 時代は $file-yusou.out
gp_TotalTransport_LatentEngy = gp_v*gp_LatentEngy
gp_TotalTransport_LatentEngy = gp_TotalTransport_LatentEngy.mean('lon').mean('sig').mean('time')
gp_TotalTransport_LatentEngy = gp_TotalTransport_LatentEngy*LatLength

gp_TotalTransport_DryStatEngy = gp_v*gp_DryStatEngy
gp_TotalTransport_DryStatEngy = gp_TotalTransport_DryStatEngy.mean('lon').mean('sig').mean('time')
gp_TotalTransport_DryStatEngy = gp_TotalTransport_DryStatEngy*LatLength


#  _  _
# [S][v] : Mean transport
#
gp_MeanTransport_LatentEngy = gp_v_tmean_xmean*gp_LatentEngy_tmean_xmean
gp_MeanTransport_LatentEngy = gp_MeanTransport_LatentEngy*LatLength

gp_MeanTransport_DryStatEngy = gp_v_tmean_xmean*gp_DryStatEngy_tmean_xmean
gp_MeanTransport_DryStatEngy = gp_MeanTransport_DryStatEngy*LatLength

# 
#                             _  _
# Vertical Mean Component of [S][v]
# gt3 時代は $file-yusou-zonal-zmean.out
gp_MeanTransportSigmaMean_LatentEngy \
  = gp_v_tmean_xmean_zmean * gp_LatentEngy_tmean_xmean_zmean
gp_MeanTransportSigmaMean_LatentEngy  \
  = gp_MeanTransportSigmaMean_LatentEngy * LatLength

gp_MeanTransportSigmaMean_DryStatEngy \
  = gp_v_tmean_xmean_zmean * gp_DryStatEngy_tmean_xmean_zmean
gp_MeanTransportSigmaMean_DryStatEngy  \
  = gp_MeanTransportSigmaMean_DryStatEngy * LatLength


#  _  _
# [S][v] z-eddy 部分
# gt3 時代は $file-yusou-zonal-zeddy.out
gp_MeanTransportZMean_LatentEngy \
= (gp_v_tmean_xmean_zeddy * gp_LatentEngy_tmean_xmean_zeddy).mean('sig')
gp_MeanTransportZMean_LatentEngy \
= gp_MeanTransportZMean_LatentEngy * LatLength

gp_MeanTransportZMean_DryStatEngy \
= (gp_v_tmean_xmean_zeddy * gp_DryStatEngy_tmean_xmean_zeddy).mean('sig')
gp_MeanTransportZMean_DryStatEngy \
= gp_MeanTransportZMean_DryStatEngy * LatLength
# これがグラフにするもの


#               _  _
# 全輸送量から [S][v] 鉛直平均部分を引いたもの
# gt3 時代は $file-yusou-totalMzonalzmean.out
# ($file-yusou.out) -  ($file-yusou-zonal-zmean.out)
#   -> $file-yusou-totalMzonalzmean.out
gp_TotalTransportM_DryStatEngy \
  =   gp_TotalTransport_DryStatEngy - gp_MeanTransportSigmaMean_DryStatEngy
gp_TotalTransportM_LatentEngy \
  =   gp_TotalTransport_LatentEngy - gp_MeanTransportSigmaMean_LatentEngy
# これがグラフにするもの.



#
# 輸送の擾乱成分 (空間変動成分)
#

# gt3 時代は $file-yusou-eddy.out だったもの
#  __ __
# [S* v*]
gp_EddyTransport_LatentEngy =  gp_v_tmean_xeddy * gp_LatentEngy_tmean_xeddy
# lon mean
gp_EddyTransport_LatentEngy = gp_EddyTransport_LatentEngy.mean('lon')
# sig mean
gp_EddyTransport_LatentEngy = (gp_EddyTransport_LatentEngy * sig_weight).sum('sig')
gp_EddyTransport_LatentEngy = gp_EddyTransport_LatentEngy * LatLength


gp_EddyTransport_DryStatEngy =  gp_v_tmean_xeddy * gp_DryStatEngy_tmean_xeddy
# lon mean
gp_EddyTransport_DryStatEngy = gp_EddyTransport_DryStatEngy.mean('lon')
# sig mean
gp_EddyTransport_DryStatEngy = (gp_EddyTransport_DryStatEngy*sig_weight).sum('sig')
gp_EddyTransport_DryStatEngy = gp_EddyTransport_DryStatEngy*LatLength
# これがグラフにするもの

# 擾乱 (時間変動成分)
# gt3 時代は $file-yusou-dist.out だったもの
# ______
# [S'v']
gp_DisturbTransport_LatentEngy = ( gp_v_tdist*gp_LatentEngy_tdist ).mean('time').mean('lon')
gp_DisturbTransport_LatentEngy = ( gp_DisturbTransport_LatentEngy*sig_weight).sum('sig')
gp_DisturbTransport_LatentEngy = gp_DisturbTransport_LatentEngy*LatLength

gp_DisturbTransport_DryStatEngy = ( gp_v_tdist*gp_DryStatEngy_tdist ).mean('time').mean('lon')
gp_DisturbTransport_DryStatEngy = ( gp_DisturbTransport_DryStatEngy*sig_weight).sum('sig')
gp_DisturbTransport_DryStatEngy = gp_DisturbTransport_DryStatEngy*LatLength
# これがグラフにするもの

#
# Drawing
#
DCL.gropn($OPT_wsn||4)
DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7)
GGraph.set_fig( 'itr'=> 1, 'viewport'=>[0.2,0.8,0.3,0.6] )


GGraph.line(gp_TotalTransportM_LatentEngy, true, "min"=>min, "max"=>max, "index"=>855, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")
GGraph.line(gp_MeanTransportZMean_LatentEngy, false, "min"=>min, "max"=>max, "index"=>298, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")
GGraph.line(gp_EddyTransport_LatentEngy, false, "min"=>min, "max"=>max, "index"=>558, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")
GGraph.line(gp_DisturbTransport_LatentEngy, false, "min"=>min, "max"=>max, "index"=>748, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")

#DCL.grcls

#__END__

GGraph.line(gp_TotalTransportM_DryStatEngy, true, "min"=>min, "max"=>max, "index"=>855, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")
GGraph.line(gp_MeanTransportZMean_DryStatEngy, false, "min"=>min, "max"=>max, "index"=>298, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")
GGraph.line(gp_EddyTransport_DryStatEngy, false, "min"=>min, "max"=>max, "index"=>558, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")
GGraph.line(gp_DisturbTransport_DryStatEngy, false, "min"=>min, "max"=>max, "index"=>748, "annotate"=>false, "exchange"=>false, "title"=>"Energy Transport")


DCL.grcls



__END__
__END__
__END__




foreach $file ('LE', SE){



#
# 作図
  $command_uchiwake = "gtcurv $file-yusou-totalMzonalzmean.out $file-yusou-zonal-zeddy.out $file-yusou-eddy.out $file-yusou-dist.out" .
     "  y=1 wsn=$wsn lay=$lay tick=F title:$file-trnsprt $options -sg:lcorner=.false.";
   if($color eq "T"){
	$command_uchiwake = $command_uchiwake .
                    " ltype=1,1,1,1 lidx=888,298,558,748 ";
   }
    if($wsn == 2){
	$command_uchiwake = $command_uchiwake .
                    " && mv dcl.ps $file-yusou-uchiwake.ps";
    }
  system("$command_uchiwake");
}