#!/usr/bin/env ruby

require "numru/gphys"
include NumRu

RPlanet  = 6371000.00000000
GasRDry  = 287
Grav     = 9.806
PI       = 3.141592653589793115997963468544185162

uvelfilename = "U.nc"
vvelfilename = "V.nc"
psfilename   = "Ps.nc"

ekefilename  = "EKE.nc"

uvel = GPhys::IO.open(uvelfilename,'U')
vvel = GPhys::IO.open(vvelfilename,'V')
ps   = GPhys::IO.open(psfilename,'Ps')

lon_weight  = GPhys::IO.open(uvelfilename,'lon_weight')
lat_weight  = GPhys::IO.open(uvelfilename,'lat_weight')
sig_weight  = GPhys::IO.open(uvelfilename,'sig_weight')

umean=uvel.average('lon')
vmean=vvel.average('lon')

eke = 0.5*( (uvel-umean)**2 + (vvel-vmean)**2 ) * ps / Grav

eke.axis(0).integ_weight=(lon_weight.val)
eke.axis(1).integ_weight=(lat_weight.val)
eke.axis(2).integ_weight=(sig_weight.val)

eke = eke.integrate('lon')
eke = eke.integrate('lat')
eke = eke.integrate('sig')

eke = eke/(4*PI)

eke.rename('EKE')
eke.set_att('long_name','eddy kinetic energy')
eke.set_att('units','J m-2')

#< Output to file >

ekefile=NetCDF.create(ekefilename)
GPhys::IO.write( ekefile, eke )
ekefile.close