# -*- coding: utf-8 -*-

#
# = エネルギー方程式の左辺を計算するスクリプト
#
# * 履歴 (新しい日付を上に書く)
#   * 2011/11/08  試行錯誤中 -> とりあえず置いときませんか?
#   * 2011/11/01  定数に単位を与えて定義する -> やったけどまだ単位がそろわない
#   * 2011/10/25  新規作成
#
# * めも
#   * 出力結果の単位が正しくないのは, 重力加速度を定数で安直に与えているから. 
#     * なんとかしたいけど方法がわからないのでひとまず放置. 
#   * 手順
#     * 各時刻での Term の領域積分の値を求める
#       * 浮力による位置エネルギーの項
#         * とある時刻, とある高度での Term を計算
#         * 高度について for を回して積分
#         * その値をさらに東西方向について積分
#       * 運動エネルギーの項と弾性エネルギーの項
#         * とある時刻での Term を計算
#         * GPhys のコマンド (.integrate) を使って領域積分
#     * 各時刻での Term の領域積分の時間微分を求める 
#     * 結果を result.dat に出力
#   # これで正しいのかはわからない
#
##########################################################################

require "numru/ggraph"
require "narray"
include NumRu

###--- 定数とか用意 ---###

dt = 10.0
tt = 20000
#tt = 40
tnmin = 0
tnmax = tt/dt
tn = 0

x1 = 0
x2 = 10000
z1 = 0
z2 = 10000

dx = 200
dz = 200
znmin = 0
znmax = z2/dz

#Grav  = 9.800000000000001    # 重力加速度 (arare5 で使っている値)
#CpDry = 1039.642343614574    # 乾燥定圧比熱 (arare5 で使っている値)
Grav  = UNumeric[ 9.800000000000001, 'm.s-2' ]        # 重力加速度 (arare5 で使っている値)
#CpDry = UNumeric[ 1039.642343614574, 'J.K-1.Kg-1' ]   # 乾燥成分の定圧比熱 (arare5 で使っている値)
CpDry = UNumeric[ 1039.642343614574, 'm2.s-2.K-1' ]   # 乾燥成分の定圧比熱 (arare5 で使っている値)

#p CpDry * Grav
#p Grav


###--- 必要な netCDF ファイルを open する ---###

gphys10 = GPhys::IO.open('../netCDF/thermal-moist_restart.nc','DensBZ')
gphys20 = GPhys::IO.open('../netCDF/thermal-moist_restart.nc','PTempBZ')
gphys30 = GPhys::IO.open('../netCDF/thermal-moist_restart.nc','VelSoundBZ')
gphys1  = GPhys::IO.open('../netCDF/thermal-moist_VelX.nc','VelX')
gphys2  = GPhys::IO.open('../netCDF/thermal-moist_VelZ.nc','VelZ')
gphys3  = GPhys::IO.open('../netCDF/thermal-moist_PTempAll.nc','PTempAll')
gphys4  = GPhys::IO.open('../netCDF/thermal-moist_ExnerAll.nc','ExnerAll')




###--- Term の計算とその領域積分 ---###

  ###--- 変数の配列を用意 ---###
  gphysIntKE    = [0]
  gphysIntElasE = [0]
  gphysIntTerm  = [0]


  ###--- 時間について for を回す ---###
  for tn in tnmin..tnmax

    time = tn * dt
    p time

    gphysDensBZ   = gphys10.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 )
    gphysPTempBZ  = gphys20.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 )
    gphysVelSoundBZ  = gphys30.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 )
    gphysVelX   = gphys1.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ).cut( 't'=>time )
    gphysVelZ   = gphys2.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ).cut( 't'=>time )
    gphysPTempAll = gphys3.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ).cut( 't'=>time )
    gphysExnerAll = gphys4.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ).cut( 't'=>time )

    ###--- 位置エネルギーの計算: Term の初期化 ---###
    gphysIntZTerm = 0

    ###--- 高度について for を回す ---###
    for zn in znmin..znmax

      height = zn * dz
      zetto  = UNumeric[ zn * dz, 'm' ]

      gphysDensBZ   = gphys10.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>height )
      gphysPTempBZ  = gphys20.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>height )
      gphysPTempAll = gphys3.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>height ).cut( 't'=>time )

      ###--- z = height での Term を計算 ---###
      gphysTerm = - gphysDensBZ * ( gphysPTempAll / gphysPTempBZ ) * Grav * zetto

      ###--- 高度について積分 (してるつもり) ---###
      gphysIntZTerm = gphysIntZTerm + gphysTerm.val * dz

    end

    ###--- 東西方向に積分 (してるつもり) ---###
    gphysIntTerm[tn] = gphysIntZTerm.sum(0) * dx

    ###--- 位置エネルギー以外の項を計算・積分 ---###

    gphysKE = gphysDensBZ * ( gphysVelX ** 2 + gphysVelZ ** 2 ) * 0.50   # 運動エネルギーを計算
    gphysElasE = 0.50 * gphysDensBZ * ( CpDry * gphysPTempBZ * gphysExnerAll / gphysVelSoundBZ ) ** 2    # 弾性エネルギーを計算

    gphysIntKE[tn] = gphysKE.integrate( 'x' ).integrate( 'z' )           # 運動エネルギーの領域積分を計算
    gphysIntElasE[tn] = gphysElasE.integrate( 'x' ).integrate( 'z' )                               # 弾性エネルギーの領域積分を計算

    ###--- テスト出力 ---###
    
  end


###--- 時間微分の計算 ---###

  ###--- 配列を用意 ---###
  gphysDtIntKE = [0]
  gphysDtIntTerm  = [0]
  gphysDtIntElasE = [0]
  gphysEnergyLeftTerm = [0]

  ###--- 時間について for を回す ---###
  for tn in tnmin..tnmax-1

    time = tn * dt

    ###--- 各時刻での Term の領域積分の時間微分の値を計算 ---###
    gphysDtIntKE[tn] = ( gphysIntKE[tn+1] - gphysIntKE[tn] ) / dt             # 運動エネルギーの領域積分の時間微分を計算
    gphysDtIntTerm[tn] = ( gphysIntTerm[tn+1] - gphysIntTerm[tn] ) / dt       # Term の領域積分の時間微分を計算
    gphysDtIntElasE[tn] = ( gphysIntElasE[tn+1] - gphysIntElasE[tn] ) / dt    # 弾性エネルギーの領域積分の時間微分を計算

    ###--- それらの和を取る ---###
    gphysEnergyLeftTerm[tn] = gphysDtIntKE[tn] + gphysDtIntTerm[tn] + gphysDtIntElasE[tn] 

  end


###--- 結果を出力してみる ---###

  ###--- 時間について for を回す ---###
  for tn in tnmin..tnmax-1

    foo = File.open("result.dat", 'a')    # 出力先のファイルを開く (a -> add)
#    foo.puts var                          # 出力 (改行あり)
    foo.print(tn*dt)                      # 時刻を出力 (改行なし)
    foo.print(",    " )                   # これがないとスペースがあかず, 前の値と次の値の境目がわからない
    foo.print(gphysEnergyLeftTerm[tn])    # エネルギー方程式の左辺を出力
    foo.print(",    " )                   # これがないとスペースがあかず, 前の値と次の値の境目がわからない
    foo.print(gphysDtIntKE[tn])    # エネルギー方程式の左辺を出力
    foo.print(",    " )                   # これがないとスペースがあかず, 前の値と次の値の境目がわからない
    foo.print(gphysDtIntTerm[tn])    # エネルギー方程式の左辺を出力
    foo.print(",    " )                   # これがないとスペースがあかず, 前の値と次の値の境目がわからない
    foo.print(gphysDtIntElasE[tn])    # エネルギー方程式の左辺を出力
    foo.print("\n")                       # 改行
    foo.close                             # 出力先のファイルを閉じる

  end