* Horizontal wave-activity flux derived by Plumb (1985, JAS)
* See (5.7) of Plumb 1985

* Used data: monthly mean data of NCEP/NCAR reanalys 1

* Monthly-mean zonal wind (uwnd), geopotential height (hgt)

* The data are available at 
* http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html

* The unit of level is [hPa]

* Monthly-mean fields in November 2014

* Basic state: zonal-mean fields 

* Perturbation: deviation from zonal-mean fields

* Level: 250 hPa


*-----
reinit

* change the directory name
* u
sdfopen /e4b/ncep/plev/monthly/nc/uwnd.mon.mean.nc

* height
sdfopen /e4b/ncep/plev/monthly/nc/hgt.mon.mean.nc

*  gas constant
define Ra=290
* earth radius
define a=6400000

define dlat = cdiff(lat,y)*3.1415/180
define dlon = cdiff(lon,x)*3.1415/180
define coslat = cos(lat*3.1415/180)
define sinlat = sin(lat*3.1415/180)
* Coriolis parameter
define f = 2*7.24/100000*sinlat
define g=9.8

* unit [hPa]
set lev 250

* set date
set time nov2014

* For drawing polar projection map
set lon -5 365

* making basic state (zonal-mean fileds)
define uclm = ave(uwnd.1,lon=0,lon=360,-b)

define zclm = ave(hgt.2,lon=0,lon=360,-b)

* deviation from zonal-mean field
define za=hgt.2-zclm

* QG stream function
define psia=g/f*za

define dpsidlon = cdiff(psia,x)/dlon
define ddpsidlonlon = cdiff(dpsidlon,x)/dlon

define dpsidlat = cdiff(psia,y)/dlat
define ddpsidlatlat = cdiff(dpsidlat,y)/dlat
define ddpsidlatlon = cdiff(dpsidlat,x)/dlon

define termx = dpsidlon*dpsidlon-psia*ddpsidlonlon
* error, found in 2013/11/25 (Thanks to Shengping HE)
*define termy = dpsidlon*dpsidlat-psia*ddpsidlonlon
define termy = dpsidlon*dpsidlat-psia*ddpsidlatlon

* "p" is normalized by 1000hPa
define coeff=(lev/1000)/(2*a*a)
*x-component
define px = coeff/(coslat)*(termx)

*y-component
define py = coeff*termy

set lon 0 360

set gxout contour
set cint 5000000
set cthick 20
set black -0.1 0.1

* QG stream-function
d maskout( psia,  abs(lat)-10)

* horizontal wave-activity flux
set arrscl 0.5 15
set arrowhead -0.3


* maskout regions where 
* (i) QG approximation is supposed to be invalid (lower latitudes)
* (ii) westerly wind speed is weak or negative
d skip(px,6,4);maskout(maskout( py , abs(lat)-10),uclm-5)

draw title (b) WAF in Nov2014 Plumb1985
printim Nov2014-P85.png