;; Example HOPS script #9
;; regrid: Regrid ocean data from velocity grid to tracer grid
;;         to compute total energy [C_p*T + 0.5*(u^2 + v^2)]

;; Open 1 year of prognostic ocean data (24 samples at half-month intervals)
;; for a 3-degree NCAR Ocean Model integration
;; (converted to CSM netCDF format using Frank Bryan's filters)

;; Mean specific heat of ocean water (J/(kg*K))
ocncp= 4000.0

;; Mean density of ocean water (kg/m**3)
ocnrho= 1000.0

hopen, "/data/large1/CSM/3x3/G03-4.24/Y0989_0990_p.nc"

;; Get surface (z=0.) T, U, V for first sample (t=1)
SST=hget("T",z=0.,t=1,subscript=1)
SSU=hget("U",z=0.,t=1,subscript=1)
SSV=hget("V",z=0.,t=1,subscript=1)

;; Convert velocities to m/s
SSU= hop(SSU, "*", 0.01, units="meters/second")
SSV= hop(SSV, "*", 0.01, units="meters/second")

;; Shift U, V to tracer grid by averaging adjacent values
SSU2= hregrid(hregrid(SSU,"x",op="zcen"), "y", op="zcen")
SSV2= hregrid(hregrid(SSV,"x",op="zcen"), "y", op="zcen")

;; NOTE: The area weights array will be deleted by the above operation

;; Compute kinetic energy on tracer grid
SSKE= hop( 0.5, "*", hop(hop(SSU2,"^",2),"+",hop(SSV2,"^",2)), name="KE" )

;; Plot kinetic energy
hplot, SSKE, fill=1

;; Compute total energy
SSTE= hop( hop(ocncp, "*", SST), "+", SSKE, name="TE", units="J/kg" )

;; Plot total energy
hplot, SSTE, fill=1

end