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

* Used data: daily 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]

* Field: 10-day averaged daily-mean fields from 11th through 20th January 1990

* Basic state ; zonal-mean fields 

* Perturbation ; deviation from zonal-mean fields

* Level; 250 hPa

*-----
reinit

* change the directory name
* u (daily-mean field)
sdfopen /e4b/ncep/plev/daily_mean/nc/uwnd.1990.nc
* height (daily-mean field)
sdfopen /e4b/ncep/plev/daily_mean/nc/hgt.1990.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

* For drawing polar projection map
set lon -5 365

* Used for judgement whether the latitude is westerly or easterly
define uclm = ave(ave(uwnd.1,lon=0,lon=360,-b),time=11jan1990,time=20jan1990)

define zclm = ave(ave(hgt.2,lon=0,lon=360,-b),time=11jan1990,time=20jan1990)

* deviation from zonal-mean field
define za=ave(hgt.2,time=11jan1990,time=20jan1990)-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 30
set black -0.1 0.1

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

* horizontal wave-activity flux
set arrscl 0.5 30

* 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 WAF in 11-20Jan1990 Plumb1985
printim daily-11-20Jan1990-P85.png