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

#
# = エネルギー方程式の浮力による位置エネルギーの項の領域積分の時間微分を計算するスクリプト
#
# * 履歴 (新しい日付を上に書く)
#   * 2011/10/20  試行錯誤 その 2, 手順を追記
#   * 2011/10/19  試行錯誤 
#   * 2011/10/18  新規作成
#
# * めも
#   * 出力結果の単位が正しくないのは, 重力加速度を定数で安直に与えているから. 
#     * なんとかしたいけど方法がわからないのでひとまず放置. 
#   * 手順
#     * 各時刻での Term の領域積分の値を求める
#       * とある時刻, とある高度での Term を計算
#       * 高度について for を回して積分
#       * その値をさらに東西方向について積分
#     * 各時刻での Term の領域積分の時間微分を求める 
#       -> 結果を result.dat に出力
#   * これで正しいのかはわからない
#
##########################################################################

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

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

dt = 10.0
tt = 20000
#tt = 2400
#tt = 200
tnmin = 0
tnmax = tt/dt
tn = 0

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

dx = 200

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

Grav = 9.800000000000001    # 重力加速度 (arare5 で使っている値)


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

gphys10 = GPhys::IO.open('../../netCDF/thermal-moist_restart.nc','DensBZ')
gphys20 = GPhys::IO.open('../../netCDF/thermal-moist_restart.nc','PTempBZ')
gphys3  = GPhys::IO.open('../../netCDF/thermal-moist_PTempAll.nc','PTempAll')




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

  ###--- 変数の配列を用意 ---###
  gphysIntTerm = [0]    # Term の領域積分値

  ###--- 時間について 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 )
    gphysPTempAll = gphys3.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
#      p height
#      p "--------------------------------------"    # ただの目印

      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 * height

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

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

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

    ###--- テスト出力 ---###
    p gphysIntTerm[tn]
    p "######################################"    # ただの目印

  end


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

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

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

    time = tn * dt

    ###--- 各時刻での Term の領域積分の時間微分の値を計算 ---###
    gphysDtIntTerm[tn] = ( gphysIntTerm[tn+1] - gphysIntTerm[tn] ) / dt        # Term の領域積分の時間微分を計算

  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(gphysDtIntTerm[tn])         # Term の領域積分の時間微分を出力
    foo.print("\n")                       # 改行
    foo.close                             # 出力先のファイルを閉じる

  end