#!/usr/bin/env ruby

require "numru/gphys"
include NumRu

#< Define constants >

GasRDry  = 287
Grav     = 9.806

#< Define filenames >

sigmadotfilename = "venus_SigDot.nc"
tempfilename     = "venus_Temp.nc"

wfilename        = "venus_W.nc"

#< Opening GPhys variables >

sigdot = GPhys::IO.open(sigmadotfilename,'SigDot')
sigm   = GPhys::IO.open(sigmadotfilename,'sigm')

temp = GPhys::IO.open(tempfilename,'Temp')
sig  = GPhys::IO.open(tempfilename,'sig')

#< Creating working variables >

imax,jmax,kmax,tmax = temp.shape

xyz_SigDot = NArray.float(imax,jmax,kmax,tmax)
xyz_W      = NArray.float(imax,jmax,kmax,tmax)

xyr_SigDot = sigdot.val
xyz_Temp   = temp.val

z_Sig    = sig.val
r_Sig    = sigm.val


#< Interpolate SigDot at sig levels >

for k in 0..kmax-1
  x = (z_Sig[k]-r_Sig[k])/(r_Sig[k+1]-r_Sig[k])
  xyz_SigDot[true,true,k,true] \
  = x * xyr_SigDot[true,true,k+1,true] + (1-x) * xyr_SigDot[true,true,k,true]
end

#< Approximate calculation of vertical wind >

for k in 0..kmax-1
  for j in 0..jmax-1
    for i in 0..imax-1
      xyz_W[i,j,k,true] \
      = xyz_SigDot[i,j,k,true]*xyz_Temp[i,j,k,true]*GasRDry/(Grav*z_Sig[k])
    end
  end
end

w = temp.copy
w.replace_val(xyz_W)
w.rename("W")
w.set_att("long_name","vertical wind")
w.set_att("units","m s-1")

#< Output to file >

wfile=NetCDF.create(wfilename)
GPhys::IO.write( wfile, w )
wfile.close