require "numru/ggraph"
require "date"
include NumRu
line = 0
qsub_int = 0.0
qsub_model = 0.0
qsub_model_exner = 0.0
cpdry = 891.0
dz =200.0
flux_top = 293.10
flux_bottom = 0.0
flux_bottom_exner = 0.0



exnerBZ = GPhys::IO.open('restart_long-ExnerBZ.nc', 'ExnerBZ')
densBZ = GPhys::IO.open('restart_long-DensBZ.nc', 'DensBZ')
ptempRad = GPhys::IO.open('../time_000000000-000080000/BS1998_PTempRad_rank000000.nc', 'PTempRad')
ptempRad = ptempRad.cut('t'=>70000..80000).mean('x').mean('y').mean('t').val
exnerBZ = exnerBZ.mean('x').mean('y').val
densBZ = densBZ.mean('x').mean('y').val
qsub = Array.new(100)

for i in 1...99 do
#  p exnerBZ[i]
#  qsub_int = qsub_int + densBZ[i] * exnerBZ[i] * cpdry * qsub[i] * dz
  qsub_model = qsub_model + densBZ[i] * cpdry * ptempRad[i] * dz
  qsub_model_exner = qsub_model + exnerBZ[i] * densBZ[i] * cpdry * ptempRad[i] * dz
end

qsub_model = dz * 0.50* cpdry *( densBZ[0] * ptempRad[0] + densBZ[99] * ptempRad[99] ) + qsub_model

qsub_model_exner = dz * 0.50* cpdry * ( exnerBZ[0] * densBZ[0] * ptempRad[0] + exnerBZ[99] * densBZ[99] * ptempRad[99] ) + qsub_model_exner

flux_bottom = flux_top - qsub_model
flux_bottom_exner = flux_top - qsub_model_exner
p "int...#{qsub_model}"
p "flux bottom ... #{flux_bottom}"
p "flux bottom_exner ... #{flux_bottom_exner}"